BNL-HET-05/10, CU-TP-1127, KANAZAWA-04-13, RBRC-486 



The Kaon 5-parameter from Quenched Domain- Wall QCD 

Y. Aoki\ T. Blum2'3, N. H. Christ^, C. Dawson^, T. Izubuchi^'^, R. D. Mawhinney^, 
J. Noaki^*, S. Ohta^-^'^, K. Orginos^ A. Soni^ and N. Yamada^'^ 

"^Physics Department, University of Wuppertal, 
Gaussstr. 20, 4^ii9, Wuppertal, Germany 
"^RIKEN BNL Research Genter, Brookhaven National Laboratory, Upton, NY 11973, USA 
'^Physics Department, University of Gonnecticut, Starrs, GT 06269-3046, USA 
'^Physics Department, Golumbia University, New York, NY 10027, USA 
^Institute of Theoretical Physics, Kanazawa University, Ishikawa 920-1192, Japan 
^Institute of Particle and Nuclear Studies, KEK, Ibaraki 305-0801, Japan 
The Graduate University for Advanced Studies (SOKENDAI), Tsukuba, Ibaraki 305-0801, Japan 

^Genter for Theoretical Physics, Massachusetts 
Institute of Technology, Gambridge, MA 02139, USA 
^Physics Department, Brookhaven National Laboratory, Upton, NY 11973, USA 

(Dated: August 10, 2005) 

We present numerical results for the kaon S-parameter, Bk, determined in 
the quenched approximation of lattice QCD. Our simulations are performed using 
domain-wall fermions and the renormalization group improved, DBW2 gauge action 
which combine to give quarks with good chiral symmetry at finite lattice spacing. Op- 
erators are renormalized non-perturbatively using the RI/MOM scheme. We study 
scaling by performing the simulation on two different lattices with a^^ = 1.982(30) 
and 2.914(54) GeV. We combine this quenched scaling study with an earlier calcu- 
lation of Bk using two flavors of dynamical, domain-wall quarks at a single lattice 
spacing to obtain Bf^ ^^^(^ = 2 GeV) = 0.563(21) (39) (30), were the first error is 
statistical, the second systematic (without quenching errors) and the third estimates 
the error due to quenching. 
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I. INTRODUCTION 

Kaon decays to two pions provided the first experimental observation of CP violation 
about four decades ago. This type of CP violation, called indirect CP violation, proceeds via 
mixing of K'^ and k'^ . Direct CP violation in kaon decays, occurring in the decay process 
itself, has been accurately measured experimentally relatively recently. Additionally, CP 
violation has now been observed in the b quark system. The standard model origin of CP 
violation is the Cabibbo-Kobayashi-Maskawa (CKM) matrix and determining the four real 
parameters that define this matrix, and looking for physical processes that are not correctly 
represented by it, is a major focus of particle physics theory and experiment. 

In determinations of the parameters of the CKM matrix, the experimental measure- 
ments of indirect (represented by the parameter e) and direct (represented by the parameter 
Re(e'/e)) CP violation in kaons should provide important constraints, but relating the ex- 
perimental values to standard model parameters requires controlled theoretical calculations. 
These calculations involve using the operator product expansion to separate the problem into 
its short-distance, and hence perturbatively calculable, components and its long-distance, 
non-perturbative parts. The short- distance effects in these processes are given by the Wilson 
coefficients of the operator product expansion [ll] and the matrix elements of the relevant 
operators in the expansion determine the long-distance parts. 

In particular, for indirect CP violation, the hadronic matrix element needed for a theo- 
retical prediction of K'^-'K^ mixing in the standard model is generally parameterized by the 
parameter Bk{h), defined by 



BM = 8-2—2 . (1) 

3J K"''K 

The four-fermion operator appearing in this expression is given 

75) d) (s 7^^(1 — 75) d) and the scale dependence, fi, enters when this operator is renormalized. 
(Here and elsewhere in this paper we will be omitting color indices for simplicity.) If one 
approximates the numerator in the definition of Bk by inserting the vacuum state to achieve 
two matrix elements of two-quark operators, one gets the value in the denominator. (This is 
known as the vacuum saturation approximation.) Since this should be a reasonable coarse 
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approximation for the Q^^'^^^^ matrix element, Bk is naturally a quantity of 0(1). 

Because Bk involves physics at low energy scales where non-perturbative QCD effects 
dominate, numerical simulations of lattice QCD provide the onlv known first principles 
method for its calculation. For recent reviews see Refs. Consequently, 
calculations of Bx have been a focus of lattice QCD simulations for two decades. Though 
achieving an accurate value for this quantity for dynamical fermion simulations is an impor- 
tant goal, a large portion of the calculations done to date are in the quenched approxima- 
tion. Early calculations used Wilson fermions, which break chiral symmetry, or staggered 
fermions, which retain a U{1) subgroup of the continuum non-singlet chiral symmetry but 
break the flavor symmetry of continuum QCD. For Wilson fermions the operator Q^^^^'^^ 
can mix with four other lattice operators, with different chiralities, at 0{a^), making precise 
calculations difficult (a is the lattice spacing) ^|^. For staggered fermions, many calcu- 
lations have been done, but the large O(a^) scaling violations of this formulation introduce 
errors in the extrapolation to the continuum limit Q]. Recently, calculations with twisted 



mass Wilson fermions 3] and improved staggered fermions |12L Il3| | have been undertaken 
to reduce these errors. 

An important improvement in the lattice techniques for calculating Bk (and other 
hadronic matrix elements) has been the development of fermion formulations which preserve 



intiy. 



and overlap fermions 



attice spacing 



15, 
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Two common 



the chiral symmetries of QCD arbitrarily well at finite 
versions of these formulations are domain-wall fermions 

ligj ]. For domain- wall fermions with their controllable breaking 



17l |. which we will use here. 



of the continuum SU^Nf)^ ® SU{Nf)ji symmetry group at finite lattice spacing, mixing of 
lattice operators, including mixing with chirally disallowed operators, is under control and 
non-perturbative renormalization techniques have been shown to work well for the relation of 
lattice operators to continuum operators [2^ Also, controlling chiral symmetry at finite 
lattice spacing removes the 0(a) scaling violations, reducing deviations from the continuum 
limit for domain-wall fermions at finite values of a. 

In this paper, we present our quenched calculation of Bk and basic low-energy hadronic 
quantities using the domain-wall fermion action and the DBW2 (Doubly Blocked Wilson in 
two-dimensional parameter space) gauge action 0, . Domain- wall fermions introduce 
a fifth dimension of length L^, with a coordinate s { < s < Lg — 1) in which the gauge 
fields are simply replicated, and produce light, left-handed quark states bound to the four- 
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dimensional boundary hypersurface (domain-wall) with s = and right-handed quark states 
on the boundary with s = Lg — 1. Four-dimensional quark fields are constructed from the 
chiral modes on the boundaries, with the residual chiral symmetry breaking controlled by the 
size of Lg. Previous works 0, 13, 13 have extensively studied the behavior of domain- 
wall fermions in quenched QCD, in particular the dependence of the residual chiral symmetry 
breaking effects on Lg. The CP-PACS Collaboration [2^ reported that the residual chiral 
symmetry breaking for domain-wall fermions in quenched QCD is markedly reduced by the 
use of a renormalization group improved gauge action (Iwasaki gauge action), which was 
also studied by the Columbia group and the smaller chiral symmetry breaking was found to 
not persist for dynamical simulations [2^. Subsequently, the RBC Collaboration found 
the further suppression of the chiral symmetry breaking by using the DBW2 gauge action, 
which was originally introduced as an approximation to the renormalization group fiow for 



lattices with a ^ ~ 2 GeV 



22 



23|. As discussed later in this paper, the Iwasaki and DBW2 



actions are closely related, differing only in the choice of a single parameter. 

Compared to the a — > value of 0.628(42) found by the JLQCD collaboration using naive 
staggered fermions [l{3], previous quenched calculations of Bk with domain- wall fermions 
have given a lower value. The CP-PACS collaboration, using the Iwasaki gauge action and 
perturbative renormalization, measured for two lattice spacings and different volumes 
and quoted a value of 0.5746(61) (191) for a ^ [s^. The RBC collaboration reported a 
value of 0.532(11) usiug the Wilsou gauge action and a smgle lattice spacing of a- = 1.92(4) 
GeV [21] and a value of 0.495(18) in a full QCD simulation with the DBW2 gauge action and 
two-fiavors of dynamical quarks at a lattice spacing of a^^ = 1.691 GeV j^ll- All of the values 
quote above refer to the B^ defined in the MS scheme using naive dimensional regularization 
and a renormalization scale fi = 2.0 GeV. This quantity has also been measured in quenched 
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33[. Here we report on calculations 



QCD using the closely related overlap formulations 
at two lattice spacings, = 1.982(30) and 2.914(54) GeV, allowing the determination of 
the a ^ value for B^- In addition, we compare with the previous RBC two-fiavor result to 
estimate the error introduced by the quenched approximation. The smaller residual chiral 
symmetry breaking for the DBW2 gauge action allows us to check whether physical results 
depend on this residual breaking. 

The remainder of this paper is organized as follows. Section |n] is devoted to a descrip- 
tion of the details of our numerical simulations and the issues of generating an ensemble of 
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configurations with the DBW2 gauge action which sample different topological sectors. Our 
results for basic quantities such as the hadron spectrum and the residual chiral symmetry 
breaking are presented in Section IIIIl We present the details of the calculation of the decay 
constants of the pseudoscalar meson in Section IIVI One of the features of our calculation 
is the use of non-perturbative renormalization (NPR) in the determination of the renor- 
malization factors needed to relate the lattice operators to their continuum counterparts. 
We deal with this topic in Section In Section IVTl we construct physical results for B^, 
compare our result with previous work and discuss the potential systematic errors for our 
result. Section gives our conclusions. Preliminary results of the calculations presented 

nn 

here have been given in Refs. |84L I85l| . 



II. DETAILS OF NUMERICAL SIMULATIONS 



Simulation Parameters 



The main results of this paper are from two ensembles of quenched configurations, one 
with = 1.982(30) GeV and the second with 2.914(54) GeV, which are generated with 
the DBW2 gauge action with f3 = 1.04 and 1.22, respectively. The DBW2 action is defined 
by 



Sg[U] 



3 



;i-8ci) J2 P[uu^u + c, J2 Rp] 



(2) 



where P[U]x,fj.u and R[U]x,fj.u represent the real part of the trace of the path ordered product 
of links around the 1x1 plaquette and 1x2 rectangle, respectively, in the /i, u plane 
at the point x and (3 = Q/g"^ with g the bare coupling constant. For the DBW2 gauge 
action, the coefficient Ci is chosen to be —1.4069, using the criteria that this action is a good 
approximation to the renormalization group flow for lattices with a^^ ^ 2 GeV 
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The 



DBW2 action is a particular choice of the class of improved actions given by plaquette and 
rectangle terms, with the Iwasaki action [s^lsil being another common choice. However, as 



was demonstrated in Ref . [2^ , the DBW2 action produces smaller residual chiral symmetry 
breaking, at a given lattice spacing and Lg, than the Iwasaki action. Our conventions for 
the domain wall fermion operator are as in Ref. |27|. 

We will also have reason to compare our DBW2 results to those obtained on quenched 
configurations generated with the Wilson gauge action, with = 1.922(43) GeV. These 
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I27I I . Table HI lists simulation parameters 



configurations were analyzed in detail in Refs. 
for the numerical calculations we present in this paper. In this table and following in 
this paper, we refer our ensembles as "DBW2 (3 = 1.22", "DBW2 (3 = 1.04" and "Wilson 
(3 = 6.0". The number of configurations used for each quantity are given, with those 
in bold denoting new calculations for this work. For the quantities where the number of 
configurations is followed by an asterisk, we employed wall-source quark propagators which 
were an average of quark propagators with periodic and anti-periodic boundary conditions 
in the time direction. These doubled the period of the correlation functions and removed 
contributions to our correlation functions from states propagating backward through the 
time boundaries of our lattices. The quark masses mj used for each calculation are listed in 
Table HH 



B. Configuration Generation with the DBW2 Action 

As mentioned previously, we have used the DBW2 gauge action for this work, since our 
earlier studies [29 1 showed this action produced a pronounced decrease in the residual chi- 
ral symmetry breaking for domain-wall fermions. This decrease occurs because the DBW2 
action gives rise to smoother gauge fields at the scale of the lattice spacing, when compared 
to other actions at the same lattice spacing. This smoothness means smaller perturbative 
contributions to the residual chiral symmetry breaking as well as many fewer lattice dislo- 
cations, where by dislocation we mean a localized change in the topology of the gauge field 
which produces eigenmodes of the domain-wall fermion operator which are undamped in the 
fifth dimension. It should be emphasized that the DBW2 action does not suppress large, 
physical topological objects, but merely small dislocations where topology is changing. This 
desired suppression of lattice dislocations has the unwanted effect of causing current heat- 
bath and overrelaxed pure gauge algorithms to sample different topological sectors quite 
slowly. 

Since topological charge changes less frequently with the DBW2 gauge action than other 
choices for the gauge action, one should check the distribution of the topological charge for 
an ensemble under study. In Ref. [2^, topology change for one of our choices of parameters, 
DBW2 with (3 = 1.04 was examined. By using many steps of a combined heat bath and 
overrelaxed algorithm, a practically acceptable frequency for topology change was observed. 
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However, here we are also interested in DBW2 lattices with (3 — 1.22. At this much weaker 
coupling and smaller lattice spacing, the time for topology change with known algorithms 
should be much longer. (We will say more about this time scale shortly.) We were forced to 
consider how to generate a distribution of DBW2 lattices at this lattice spacing, including 
configurations of different topology. The different topological sectors for DBW2 should be 
distinguished by large, physical topological objects. 

Consider the partition function for quenched QCD in a particular topological sector with 
topology n. For the DBW2 action, this is exphcitly given by 

Z^L2 = / MC/]e-^— (3) 

Jn 

where only gauge fields with topology n enter. Of course this requires a precise definition 
of the topology of the gauge field and one could use the domain-wall fermion operator, 
for arbitrarily large L^, to determine this. With the assumption that current algorithms 
do not change the topology for DBW2 lattices with /3 = 1.22, a thermahzed lattice with 
topology n will remain in this topological sector. This makes it straightforward to measure 
the expectation value {0)n of an observable O in the topological sector n from a starting 
lattice with a given topology. Since Z = J2n Zdbw2-> "^^ have 

l^m ^DBW2 

and we require the ratio '^dbw2/'^dbw2' configurations in topological sectors n and m, 
to determine (O). Unfortunately, this ratio is not simple to determine. 

We can, however, approximate this ratio using a corresponding ratio with the standard 
Wilson gauge action at the same lattice spacing and volume, i. e. 

y{n) y{n) 

^DBW2 ^ ^Wilson /r\ 
y{m) ~ y{fh) V / 

^DBW2 ^Wilson 

where the n and fh variables denote that topology in the Wilson case is determined only 
from the long-distance features of the gauge fields, n = n and m = ffi. This approximation 
is motivated by the expectation that large-scale, physical topological fluctuations will be the 
same for any two theories that differ only in their cut-off behavior at short distances. This 
is a standard statement of field theory. An alternative way of expressing this is to consider 
the topological susceptibihty. When averaged over all topological sectors, it should give the 
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same result for Wilson and DBW2 actions at weak coupling, provided the susceptibility is 
renormalized correctly in both cases. Renormalization improvement removes the dependence 
on the ultraviolet parts of the theory and is equivalent to using only large-scale topological 
objects in determining n and fh. 

To use Eq. inito generate an ensemble of DBW2 configurations with the correct large-scale 
topological fluctuations of quenched QCD for DBW2 with (3 = 1.22, we proceed as sketched 
in Fig. H We first generate gauge configurations using the standard Wilson gauge action 
with P = 6.25, denoted by proceeding down the vertical dashed line in Fig. ^ By using 
nip = 770 MeV as input, these lattices have a lattice spacing of 2.808(49) GeV with = 8 
and Ms = 1.7 for 50 configurations. Alternatively we can use the heavy quark potential 
to set the scale. With the parameterization in Ref. [3| and a value for the Sommer scale, 
tq, of 0.5 fm, we obtain = 3.12 GeV for the lattices with Wilson gauge action, which 
is close to the lattice s pac ing of = 3.09(2) GeV as determined by the heavy quark 
potential for f3 = 1.22 39]. Every 10,000 heatbath sweeps (the white circles in Fig. P), 
a Wilson lattice is saved to use as a starting point for a DBW2 evolution (the horizontal 
lines in Fig. Q). Assuming the initial Wilson gauge configurations effectively sample the 
different topological sectors, we are then beginning each DBW2 evolution from a starting 
configuration which reflects the appropriate large-scale topological distribution of quenched 
QCD. The parameters of the DBW2 evolution are chosen to yield the same lattice spacing 
as in the Wilson evolution. We evolve using DBW2 to reach an equilibrated ensemble for 
DBW2, but in a particular topological sector, assuming that the topology does not change 
during the DBW2 evolutions. We can then average observables from the different DBW2 
evolutions together, since the probability of each topology appearing is controlled by the 
Wilson action and Eq. El relates the Wilson and DBW2 quantities. 

For observables that are not sensitive to short- distance topological features, the algorithm 
above should provide a good approximation to a full DBW2 ensemble average. The appro- 
priately renormalized topological susceptibility is an example. However, an uncontrolled 
approximation in this algorithm is how the topology in the initial Wilson action changes 
as the DBW2 evolution thermalizes. Small size topological defects in the Wilson lattice 
are suppressed by the DBW2 action. Defects which are removed, leaving the topology un- 
changed, cause no uncertainty. Defects which are removed by becoming larger in size can 
contribute to the topology as determined on large distance scales and make a topology in 
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the DBW2 case appear with an incorrect probabihty. 

Overall, we used 53 initial Wilson configurations (the open circles in Fig. ^ and 53 
subsequent independent DBW2 evolutions. To check the distribution of topological charge 
Qtop for our ensembles, we have calculated it for co nfigu rations indicated by the open and 
filled circles in Fig. ^ a la the MILC Collaboration |3ll. l4fll|^. Figure El shows the result in 
the form of a time history and distribution of the topological charge for three ensembles of 
53 lattices: one with no DBW2 sweeps (the initial 53 Wilson configurations); the second 
generated from the initial Wilson configurations with 5,000 DBW2 over-relaxed/heatbath 
sweeps; the third generated with 10,000 DBW2 sweeps. These appear in Fig. Elfrom the top 
to bottom panels, respectively. It is apparent from the figure that the topology, as measured 
using the smearing technique of Ref. j40], changes little during a DBW2 evolution. During 
the 10,000 sweeps done for these 53 different evolutions, 21 evolutions changed topology 
once and the others did not. Of the 21 evolutions which changed topology, 19 changed 
in the first 5,000 sweeps. Averaging over each ensemble, we obtain (Qtop) = 0.18(43) (0 
sweeps), 0.36(40) (5,000 sweeps) and 0.40(40) (10,000 sweeps). Taking the last two ensembles 
together, our DBW2 ensemble has (Qtop) = 0.38(29). 

For the physical observables reported in this paper (the DBW2 (3 = 1.22 column of 
Table |l}, we collected one or two configurations from either 5,000 or 10,000 sweeps in each 
of the 53 DBW2 evolutions. To check for thermalization effects in the DBW2 evolutions, 
we measured the residual quark mass, whose definition is given in Eq. |Hl and pseudoscalar 
meson mass, mps, after every 1,000 over-relaxed/heatbath sweeps with = 8 and M5 = 1.7. 
The results for rrij-es from the first 20 configurations in the direction of the Wilson sweep 
are shown in Fig. El for different values of the quark mass, fnf = 0.02, 0.03 and 0.04 (open 
symbols) and those in the chiral limit from the linear fit (filled symbols). The vertical axis 
is divided into two parts, since the residual mass for the initial (Wilson) lattices is about an 
order of magnitude greater than for the DBW2 lattices. No thermalization effects are visible 
after the first 1,000 sweeps, but clearly the result for rrires shows that lattice dislocations 
have been markedly reduced. Figure IH shows the same kind of plot for mps. Again, no 
thermalization effects are visible here, after the first 1,000 sweeps. From this we see that 
DBW2 lattices separated by 5,000 sweeps are thermalized. 



We thank the MILC CoUaboration for their code which was used to compute the topological charge. 
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We thus conclude that we have generated a thermahzed distribution of DBW2 lattices, 
with a distribution of large-scale topological features that is a good approximation to the 
exact distribution. Our strategy is based on the assumption that a modification of the ultra- 
violet properties of the theory, such as the RG- improvement of the action, does not change 
the infrared properties of the theory, such as topology. The measurements discussed here 
show that the approximation is working quite well. We believe this approach to be superior 
to either working in a single topological sector, since our lattice volumes are not large, or 
averaging randomly over topologies, since that ignores the underlying QCD dynamics. Of 
course, having a pure-gauge updating algorithm which rapidly samples different topological 
sectors at weak coupling would be a better solution, but with this current approach we now 
turn to the measured meson masses. 



In this section we discuss the calculation of meson masses and use them to determine 
values for the basic parameters of the simulation, such as the residual quark mass, m^es, the 
lattice spacing, a, and the bare quark mass which produces a pseudoscalar meson with the 
physical kaon mass when made from two degenerate quarks. 

We have measured the wall-point and wall-wall two-point correlation functions 



where (prit) and Xr(^) are quark bilinear interpolating fields with the Dirac spinor structure 
r = V^,A^,S and P. The quantity (f)r{t) is a local, bilinear field, summed over a spatial 
volume at fixed time t, that is 4>r(t) = J^x'li'^jt) rq'(x, t). In contrast, Xr(t) is a spatially 
non-local bilinear defined on a three-dimensional volume at fixed time. In particular, Xr{t) = 
J2x,y ^(x, t) r g(y, t) for Coulomb gauge fixed quark fields q on time-slice t. 

As explained in Section III Al we first compute quark propagators with both periodic 
and anti-periodic boundary conditions in the time direction in evaluating the correlation 
functions in Eqs. El and [7| from certain ensembles. To extract the masses and ampli- 
tudes for states entering these correlation functions, they are fit to the hyperbolic functions 
Asmh{m(t — to — T/2)) or Acosh{m(t — to — T/2)), depending on the symmetry of the 



III. VECTOR AND PSEUDOSCALAR MESON MASSES 




(6) 
(7) 
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interpolating fields being used. Here A and m are fitting parameters and T is twice the 
time extent of the lattice, due to our choice of quark propagators. We note that in our 



earlier quenched work [22, |22( for the Wilson /5 = 6 and DBW2 /3 = 1.04 actions, we did 
not calculate two-point correlation functions on doubled lattices. A further improvement 
in the present work is the use of two-point correlators which are the average of two point 
functions obtained from each of the two sources that were introduced to compute the three 
point functions. The data obtained in this way are marked with *'s in Table H] 

For domain-wall fermion simulations the finite extent of the fifth dimension produces 
chiral symmetry breaking effects in the low-energy QCD physics represented by the domain- 
wall fermion modes localized on the boundaries of the fifth dimension. We first turn to a 
determination of this residual chiral symmetry breaking since it is of both intrinsic interest 
and is needed for all extrapolations to the zero (renormalized) quark mass limit. As has 



been discussed extensively in Ref. 



27l |. for low-energy QCD the residual chiral symmetry 



breaking will appear as a small additive quark mass, denoted by mres, that represents this 
symmetry breaking in an effective-field theory formulation of QCD. 
We calculate mrps from the ratio of correlators 





J5,(x,()fJ.(0) 




(0 


MmUo) 





rrires = r 1 ^ — '- , (8) 

where ^ = for Wilson p = 6.0 and DBW2 p = 1.04 and ^ = x for DBW2 /? = 1.22 
depending on an unessential convenience in our numerical simulation. J^g is a pseudoscalar 
density located at the mid-point of the fifth dimension jl?!. |27|. In Fig. EJwe show the ratio 
on the right-hand side of Eq. IHlas a function of t for DBW2 f3 = 1.22. From the bottom to 
top, the panels in Fig. Elshow data in order from the lightest mj to heaviest. In choosing 
the lower bound of the fitting range tmin, we only need a time separation large enough to 
remove any contribution from the unphysical, off-shell, five-dimensional states in the DWF 
theory. We do not need to suppress legitimate excited states of QCD whose contributions 
will also obey Eq. |H1 As suggested by the figure, we could chose tmin as small as 4-5 lattice 
units. For somewhat arbitrary reasons we used the value tmin = 14 which has the effect of 
somewhat increasing the resulting error on mj-es- 

Since the finite Lg effects represented by rrires arise from short distances, mres should 
depend only weakly on the quark mass and is usually evaluated in the limit 0. A 

plot of the values of function of m/ is shown in Fig. IHl Included in the figure is 
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a linear fit to the data, which we extrapolate to rrif = to find mres = 0.9722(27) x 10~^. 
The extrapolation changes the value of rrij-es by less than 1%. The coefficients of the linear 
fit are given in the first row of Table IIIIl The second and third rows give results for mres 
previously calculated for other ensembles that will be used in this paper. We note that mres 
for the Wilson gauge action with jS = 6.0 was obtained from simulations with a single quark 
mass of TUf = 0.02. 

In a previous paper on simulations with domain- wall fermions j27|, we included an 
extensive investigation of the infrared pathologies that occur in the — limit of quenched 
domain-wall fermions, or any other fermionic formulation that preserves the continuum 
global symmetries of QCD. In this paper, the focus is primarily on physics at the kaon scale, 
but we will briefiy check that our earlier observations are consistent with the data at this 
weaker coupling. To this end, we now detail our methods for determining physical values 
from pseudoscalar correlators. 

In extracting the low-lying masses and amplitudes from the correlators C^(t, to) and 

topological near-zero modes. We de- 



Cpp (t, to) it is important to minimize the effects o 



scribed a number of ways to approach this in Ref. 22| and noted that the effects of topo- 
logical near-zero modes decrease as the source-sink separation is made larger. Here we 
couple this observation with our measurement of correlators on doubled lattices to minimize 
the effects of topological near-zero modes by choosing a relatively large value of tmin when 
extracting masses and amphtudes. 

We now describe our results for the pseudoscalar masses and leave the detailed discussion 
of the pseudoscalar decay constants to the next section. In Fig. Owe show the effective mass 
of the pseudoscalar meson mass as a function of time for the DBW2 f3 = 1.22 data set 
obtained from the point- wall correlator C^^. A fine plateau showing no apparent excited 
state contamination extends from t = 10 to t = 38. Note the source is located at t = and 
the symmetric mid-point of our time-doubled lattice is t = 47.5. In this paper we use an 
analytic formula to determine the effective mass, rrics, from three time separations: 



meff = ln(r(t) + ^r(t)2 - 1) (9) 

where the ratio r{t) is given by: 
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and C(t) represents one of the two-point correlators defined in Eqs. iniand[7| 

Since tlie decay constants will be calculated and compared using both C^^^ and C^^ m 
the following section, we want to compare the amplitudes determined from fitting the two 
different correlators, while keeping the pseudoscalar mass the same. To achieve this, we 
extract a common value of mps from these correlation functions through a simultaneous fit 
which minimizes the given by 



X 



E 



C^^^it, 0) - A^i^ sinh(mps(t - T/2)) 



+ 



C^^(t, 0) - AZ cosh(mps(t - T/2)) 



PP 



(11) 



where a(t) is the iackknife error of the correlator at t. 

n 

In Ref. j27|], it was found that, for Wilson (3 = 6.0, the zero mode effects seen by comparing 
scalar and pseudoscalar correlators become small for source-sink separations of 10 lattice 
spacings. To mitigate these effects in the analysis of this work, we chose the fitting range 
(tmin, tmax) to be (12, 19) for DBW2 f3 = 1.04 and Wilson /3 = 6.0, both of which correspond 
to ^ 2 GeV. For DBW2 (3 = 1.22 (a^i ^ 3 GeV), we use (18, 31) so that tmin = 12 and 
18 corresponds to a similar distance in physical units for both lattice spacings. Results for 
mps with degenerate and non-degenerate masses for DBW2 /5 = 1.22 are given in the fourth 
column of Table HVl Results of the same analysis on the doubled lattice for DBW2 f3 = 1.04 
and Wilson f3 = 6.0 are listed in the sixth column of Table IVl This table contains previous 
values for mps from the point-point correlator C^^ in the fifth column for comparison. The 
central values agree within the quoted errors, with smaller errors for the results computed 
on larger ensembles. For the DBW2 action with /3 = 1.22, we have also calculated mps 
from Cp^^, C^^ and C^^ separately and the results are consistent with those from the 
simultaneous fit obtained by minimizing the given in Eq. ^2 within the quoted statistical 
error. 

We now investigate the chiral limit. In Ref. j^] 

{mf + TJly-es) 



mpg is fit to the form 



mpg = a.„{mf + m^ 



6 In 



(12) 



with the cutoff = 1 GeV. This expansion includes the pathologies of the quenched 
approximation 5m f Inm/ making it more difficult to extrapolate to the chiral limit whereas 
only weaker m^jlnmf terms appear in the full theory 4l|, |42, |43|. The coefficient 6, as 
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well as the other fitting parameters, must be obtained from the data itself. Excluding the 
next-to-leading order (NLO) coefficient (6 = 0), the fitting parameters corresponding to four 
varieties of range are given in Table IVTl for DBW2 j3 = 1.22. From this table, one sees 
consistent small values for 6 with large uncertainties. Similar results appear in Ref. 2£|] for 
DBW2 P = 1.04 and in Ref. for Wilson /3 = 6.0, though in the latter a slightly different 
parameterization was used. To get meaningful result for b ^ 0, one must include the full 
covariance matrix which leads to rather poor values of for the fits given in the last two 
rows of Table IVII The third and fourth rows of the Table shows the result of excluding the 
lightest mass point to avoid possible near-zero mode contamination which is present in the 
quenched approximation 

Since for small quark masses, the quenched chiral logarithm dominates the term quadratic 
in quark mass, it may be most reliable to determine delta by using only light quark masses 
and setting 6 = 0. Otherwise, the non-linearities due to the quenched chiral logarithm are 
being offset by the quadratic term and this leads to a marked change in the value for a-,,, as 
can be seen from Table IVTl Of course, even smaller quark masses are required to completely 
justify the omission of the quadratic term. 

Topological near-zero modes of the Dirac operator, which are not suppressed in the 
quenched approximation, can give rise to a non-zero value of mps in the chiral limit |27|. 
This effect manifests itself as a finite volume effect since the density of such modes decreases 
as ~ 1/VV. Since we focus on the region around ms/2 for the determination of the kaon 
S-parameter, in what follows we ignore these effects which only become important near 
the chiral limit and use a simple linear fitting function for mpg. A definitive study of 
the parameter 6, which is not the goal of this work, requires smaller quark masses, larger 
statistics, and larger physical volumes than have been used here. 

We determine the lattice spacing from the vector meson mass. To improve the qual- 
ity of the signal, the vector correlation function is averaged over all spatial polarizations, 
ELiC^^VS- We set the lattice spacing by extrapolating my in lattice units to the chiral 
limit {jnf = —m^es) and equate this value to a x 770 MeV, the same procedure we followed 



in our previous quenched studies 
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27 



Figure IHl shows the vector meson effective masses with degenerate quark masses for DBW2 
(3 = 1.22. The lines denote results of fits to the hyperbolic function mentioned above and 



indicate the fitting range, central value, and magnitude of the jackknife error for each. 
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Results for mv for each value of rrif are collected in the second column of Table HVl The 
values of the ratio mps/my are also included in the third column of this table. Similar 

n n 

analyses were carried out in Refs. |29l] and [221 DBW2 (3 = 1.04 and Wilson p = 6.0, 
respectively. Results for niv and rrips/fnv are quoted from these references in the second 
and third columns of Table IVl 

Examples of the chiral extrapolation for my and mpg with DBW2 (3 = 1.22 are shown 
in Figs, ini and respectively. These figures contain masses obtained with non-degenerate 
quarks masses mi and m2 as well. Data are plotted as rrif = (mi + m2 + 2mres)/2. They 
appear to lie on a smooth line joining the degenerate mass points. We take linear fitting 
functions for the pseudoscalar mass-squared and the vector mass, 

mpg = Co + cf (m/ + mres), (13) 
my = Cq + c^{mf + rrij-es). (14) 



Values of these parameters are tabulated in Table UTTl These parameters were determined by 
minimizing an expression for in which off-diagonal terms in the covariance matrix were 
omitted. Such an uncorrelated fit will give an valid result but, if substantial correlations 
are present, will result in a somewhat larger statistical error and an unusually small value 
for x^- Singularities in the full covariance matrix prevented our including the off-diagonal 
terms. Note, here mres is not a free parameter, but is fixed to be the central values of the 
results from Eq. |H1 Because of the small values of m^es in our simulations, we do not take 
the non-zero value of Cq as an indication of explicit chiral symmetry breaking. As explained 
above and in detail in Ref. ^31, both unsuppressed topological near-zero modes of the Dirac 
operator in quenched simulations and neglecting the quenched chiral logarithm can make 
demonstrating that m^ is zero in the chiral limit difficult. 

We have also computed the values of the bare strange quark mass m^ and J-parameter 
for DBW2 ensembles and results are summarized in Table I VI II The strange quark mass is 
found from 

y^c^ + cf(m,/2 + mres) _ rriK _ 495 MeV 

+ cY{mud + m,es) ~ 770 MeV ^ ^ 

Here we simply set rriud = — '^res- The presence of a non-zero intercept Cq for mpg in Eq. fT!^ 
makes a more precise determination of rriud difficult. The overall effect of this approxima- 
tion is an ^ 1% error in the determination of rrip. We extract the kaon i?-parameter by 
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interpolation to the physical point, = 495 MeV, which is equivalent to evaluating our 
data at the point mj = 171^/2 determined from Eq. ^1 The J-parameter is defined as in 
Ref. Q 

dmv cX , . 

J = -, rriK* = -pTTiK*, (16) 

amps cf 

where mx* = Cq + c({ms/2 + mres)- 

The chiral symmetry of domain- wall fermions suppresses 0{a) discretization errors in 
low-energy observables, so the leading error is expected to be O(a^). Thus, to study the 
scaling dependence of the J-parameter, we plot our results and the previous one for DBW2 
13 = 0.87 Q in Fig. ^2 as a function of [GeV ^] . These results are consistent with 
each other, showing that discretization errors are indeed small for this quantity. However, 
the quenched J-parameter value is much smaller than the experimental value, mx* {rriK* — 
^p)/i^K - ^l) = 0-48, by about 30%. 

In making the comparison described above and shown in Fig. we should emphasize 
that in addition to the quantified variation in lattice spacing, these points also correspond 
to varying gauge actions (Wilson and DBW2) and to different values of the domain-wall 
height, M5. Since the coefficient of the O(a^) correction to the J parameter can depend on 
both the action and M5, we should not attempt to fit the points in Fig. to a single linear 
term in a?. However, the agreement between these various values of J certainly suggests 
considerable independence of the lattice spacing over a large range of lattice scales. 

Before concluding this section, we present values for m|g in physical units [GeV^] in 
Table I VI HI In the following sections, these values will provide a physical horizontal scale, 
when quantities such as the decay constants, K-K matrix elements and Bk are plotted as 
functions of quark mass. When mpg is used in this way, the chiral limit will be identified as 
the point m|g = 0. However, as was mentioned above and will be discussed further later, 
the slight difference between the points m|>g = and rrif = —mres when a simple linear fit is 
done (likely arising from the effects of near zero modes and neglecting the quenched chiral 
logarithm term) will introduce systematic errors at or below the 1% level. 
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IV. PSEUDOSCALAR MESON DECAY CONSTANTS 

In this section we discuss the calculation of the pseudoscalar decay constant /ps- In our 
Euclidean conventions, the definition of this quantity is 

(OIA4IPS) =/ps-mps, (17) 

where |PS) stands for the zero-momentum pseudoscalar state. The experimental values for 
the pion and kaon states are ~ 130 MeV and fx ~ 160 MeV. 

We consider three separate lattice transcriptions of Eq. El In each case, rather than 
use the conserved axial current, we make use of the local, fiavor non-singlet axial current. 
These two quantities are related by a finite renormalization constant, and we first discuss 
the extraction of the "bare" value of /ps, to which this factor has yet to be applied. The 
three formulae we use to extract the decay constant are 

AAiP 

ffi = — ^i^=, (18) 

2 -^ww 



'PS 



f(2) 

/.(3) _ / -^>^w ^ "^pw 

/ps — W T7~" ^ 



(19) 
(20) 



t>ti 



where C^^^ and are point-wall and wall-wall, axial-pseudoscalar and pseudoscalar- 
pseudoscalar correlation functions discussed in the previous section, mpg is the common 
mass extracted from these correlation functions, and 



Kf = (0 \M PS) (ps xi 0) , (21) 



2mps 

<.^-^l(PS|XP|0)r, (22) 

are expressed as the product of correlation function amplitudes extracted from the same fit. 

We take /pg* for our final value, using the same fitting range used to extract the pseu- 
doscalar mass in the previous section. However, by also calculating fp^ and /pg* we are able 
to study the size of the systematic error that may arise from our choice of fit which may 
affect, for example, the amount of excited state contamination in the result. Another source 
of systematic error is unsuppressed topological near- zero modes of the Dirac operator in the 
presence of quenched gauge configurations. This contamination is elucidated in detail in 
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Ref . [^^I , the salient points being that the contamination depends on the time separation of 
the operators, the particular construction of these operators (for example, whether a wall or 
point source is used), and the value of quark mass, lighter quarks showing a larger effect. The 
comparison between /pg* and /pg* is particularly interesting because the expression we use to 
calculate Bk essentially contains /pg in the denominator, /pg tends to be contaminated by 
the topological near-zero mode since it is extracted from time slices closer than the case of 
/pg . As such, we use the same source positions and fitting ranges for the extraction of /pg* 
as we do for Bk, (^1,^2) = (7,41) for DBW2 p = 1.22 and (5,27) for p = 1.04 and Wilson 

P = 6.0, with final values given by the fit over the ranges 19 < t < 29 and 14 < t < 17, 

(3) 

respectively. For the purpose of this comparison we choose the fitting range for /pg to be 
the same as that chosen for /pg*. 

Table IIXI shows the results for all three values of the decay constant for each quark 
mass. As can be seen, all these quantities agree within their quoted statistical errors. 
However, the statistical fluctuations of these quantities are correlated, so differences between 
these quantities may be resolved to a much higher precision than the quantities themselves. 
Figure El shows (/pg* — fps)/fps ^ function of m|g for each ensemble. The central value 
of this difference monotonically increases and becomes roughly one percent at the kaon mass 
point mpg = rriK = 495 MeV. Although data except at the lightest masses are consistent 
with zero within one standard deviation, it is important to note that the same behavior 
is observed for all independent ensembles. While it may suggest the effects of topological 
near-zero modes, more statistics are required for further study. In Section IVIB| we will 
discuss systematic error of Bk which may originate from the ambiguity of /ps stemming 
from the different methods of extraction. 

As mentioned previously, in the rest of this paper we employ /pg* as our result for the 
decay constant. In Figs. Eland El effective values of /pg* are plotted versus time for DBW2 
P = 1.22 and 1.04, respectively. While no significant time dependence is seen within the 
choice of fitting range for (3 = 1.04, we find a monotonic decrease for smaller t for (3 = 1.22. 
However, the width of this variation is within the statistical error. 

The renormalization constant of the local, flavor non-singlet, axial current is determined 
following the method in Ref. j27|. As an example, the results for each quark mass for 
DBW2 (3 = 1.22 are shown in Fig. El is defined in the chiral limit which in practice is 
determined from a linear extrapolation to m/ = —rrires- Results for Za are summarized in 
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Table H 

To determine the physical decay constants f^^ and fx we make use of NLO quenched 
chiral perturbation theory which suggests a simple linear quark mass dependence in the 



case of degenerate quarks |4^ 3]. The results of these linear fits of the bare value /ps are 
listed in Table 1x1 In the same table, the lattice values of decay constant /^^^**^ and /j^*^**^ are 
listed. They are obtained by the extrapolation to the point mpg = = 135 MeV and the 
interpolation to mps = mx = 495 MeV, respectively. In particular for the DBW2 results, 
we note that the lattice scale dependence of the renormalized value Z^fps comes mainly 
from the scale dependence of Z^. In Fig. UHlthe renormalized decay constants Z^-fps are 
shown along with linear fits for DBW2 (3 = 1.22 (solid line) and (3 = 1.04 (dashed line). One 
observes agreement between DBW2 /? = 1.22 and Wilson P = 6.0 values, but a discrepancy 
of roughly two standard deviations between DBW2 (3 = 1.22 and (3 = 1.04. 

Renormalized decay constants /^r = ZaJ^^^^^ and fx = ZaJk^^^^ and their ratio fxl f-K = 
yjiatt) ^y.^iatt) g^j^g listed in Table IXTl Since the last quantity contains neither the statistical 
errors of nor the scaling dependence of Za-, it may allow a more accurate study of scaling 
than the individual decay constants. Putting together our /k/Ztt results with previous ones 
for DBW2 P = 1.04 and 0.87 Q, 

which were obtained from a point-point correlation 
function on a non-doubled lattice, we find the scaling shown in Fig. El At f3 = 1.04, our 
result is within the statistical error of the previous one. We extrapolate our DBW2 results 
and the previous (3 = 0.87 value (filled symbols in the figure) linearly in a^, (see the discussion 
at the end of the last section) and find the continuum limit value fx/ fix = 1.098(13). A 
constant fit yields a consistent value with x^/dof = 0.91. These results are roughly consistent 
with the one-loop analytic result ~ 1.07 of quenched chiral perturbation theory 43[ and differ 
from the experimental results by several standard deviations. Similar plots are obtained for 
the individual decay constants as shown in Fig.^Jand results of the continuum extrapolation 
both with a linear and a constant fit are listed in Table IXII We again find consistent 
continuum values within the errors for both types of extrapolation. 

V. NON-PERTURBATIVE RENORMALIZATION FOR Bk 

The value for Bk depends on both renormalization scheme and scale. In this work, 
we will be quoting our final answer renormalized in the MS scheme at 2 GeV using the 
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non-perturbative renormalization (NPR) technique of the Rome-Southampton group 
This technique has been found to be very successful, particularly when used in conjunction 
with domain-wall fermions, in which context it has been applied^ the renormalization of 



fermion bilinear operators (20, |4a|, Bk in previous calculations j2l|, |4a, |42|, and the four- 
quark operators in the AS* = 1 effective Hamiltonian j2l|. 

In the continuum, the parity-even operator of interest for the calculation of B^, 

Ovv+AA = {s-ff,d){s-ff,d) + is-f5-ff,d){s-f5-ff,d) , (23) 

renormalizes multiplicatively. However, in a regularization in which chiral symmetry is 
broken, this operator may mix with four other four-quark operators; we must then solve the 
renormalization problem on the following basis of five operators: 

Ovv±AA = (s7^rf)(s7^c/) ± (s757^rf)(s757^d), (24) 
Oss±pp = isd){sd)±{s-f5d){s-f5d), (25) 
Ott = {scr^^d){sa^^d) . (26) 

When using domain-wall fermions such mixing with wrong chirality operators should be 
strongly suppressed. However, as a consequence of their different chiral structure, chiral 
perturbation theory predicts that the contribution of these operators to Bps will diverge in 
the chiral limit. As will be discussed in Section FVI Bl even for the (relatively large) masses 
at which we are working, the contributions to Bps from these wrong chirality operators are a 
few dozens of times larger than the one from Ow+aa, and so even a relatively small mixing 
coefficient may become numerically important in the final answer. In this work we address 
this problem by presenting both a theoretical argument to estimate the size of such terms, 
and the results of a numerical study. The question of possibly large mixing with wrong 
chirality operators in a domain wall fermion calculation of Bj^- was raised in Ref. 0. The 
discussions presented here are intended to resolve this issue. 

We may theoretically estimate the size of the mixing coefficients between Ow+aa and 
the wrong chirality operators by applying the spurion field technique introduced in Ref. j^l ■ 
The details of this analysis are presented in Appendix El here we merely quote the result 
that, if we consider only the effects of the explicit chiral symmetry breaking of domain 
wall fermions, these mixing coefficients occur at 0(m^gj. This represents a suppression 
factor which, in the presented calculations, is in the range 10^^ - 10^®, well below the level 
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which could make a significant contribution to our result. If we also take into account the 
fact that we are working at finite quark mass, we would expect the leading contribution 
to the mixing coefficients to be of 0(mj). Such contributions will appear in any lattice 
formulation, and should be the dominant contribution from wrong chirality operators in 
any realistic calculation of Bk using domain-wall fermions. In the following we numerically 
estimate the size of such terms. 

We calculate the renormalization coefficients in a two step process: first we calculate the 
renormalization factor on the lattice in the RI/MOM-scheme, then make use of a continuum 
running/matching calculation to convert this into the MS-scheme. In this way, we avoid 
the use of lattice perturbation theory for which the convergence properties are problematic 
(requiring the use of mean field improvement). To apply the NPR technique we work in 
Landau gauge and construct the amputated n-point correlation functions of the operators 
of interest, p^^)^^**^ with external quark lines carrying large, off-shell momenta. The renor- 
malization conditions may then be applied by requiring that suitable spin-color projections 
of this correlation function are equal to their tree case value. To be precise, for a general 
mixing problem involving m operators, Oi {i G 1, ■ ■ ■ ,^12), each made up of n quarks, we 
define the renormalized operators by the equation 

= * , (27) 

construct the amputated n-point correlation function, and apply, at a fixed configuration of 
external quark momenta, the condition 

where represents the application of a particular spin-color projection and a subsequent 
trace. (For details see Appendix^) For a mixing problem with m operators, m independent 
projection operators need to be applied, however the precise choice will not effect the final 
renormalization factors. This is in contrast to the particular choice of gauge, quark mass, 
and the configuration of external quark momenta: these must be specified to fully define 
the renormalization condition (in all cases, we define our renormalization conditions to be 
in the chiral limit). 

To calculate a renormalized value of B^, the ratio of renormalization factors Zq(a.s=2)/Z\, 
which we refer to as Zbj^, is required. As such, we calculated the amputated 3-point and 5- 
point Green functions for O = and O = Q'^^^^'^\ in each case using the same magnitude 



T(n)tree 



(28) 
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of momenta for all external quarks. As will be explained below, this is not the optimal choice 
of momenta for the lattice calculation. However, it is the only momenta configurations for 
which perturbative calculations of the matching between the RI/MOM- and MS-schemes 
exist. For the axial- vector operator the spin-color projection is achieved by tracing with 
757^. The ratio of renormalization factors, Z~^Zai is then obtained from 

Z,'ZM^^^'"\p: «-^)757.] = -48, (29) 

which is valid to O(a^) even with non-zero mj-es The details of the projection operators 
used and derivation of the counterpart of Eq. 1221 for the AS* = 2 operators are given in 
Appendix El in the following we will simply refer to the projected, amputated n-point 
correlation functions in terms of the matrices A and A^, defined as: 



A- = P- 



platt 



AT. = P 



ptree 



(30) 
(31) 

Using these definitions Eq. |2H1 takes the form of a matrix equation 

Z-2z,,A,fc = AT^fc . (32) 

In the first step of the numerical calculation, we calculate the quark propagator in co- 
ordinate space for each m f listed in Table |nl After performing a Fourier transformation 
into momentum space, we produce the quark propagator G{piatt), the treatment of which 
is described in Appendix iBl We use a set of integer momenta {nx,ny,nz,nt) on the lattice 
defined as 

(27r 27r 2n 2n \ 
Piatt = -j-n^, l-^y^ -f-nz, -f-rit , (33) 
\Lx Ly Lz j 

where = Ly = Lz = 24 and Lt = 48 for DBW2 f3 = 1.22, and L^ = Ly = Lz = 16 
and Lt = 32 for DBW2 (3 = 1.04. We used 448 combinations of {nx,ny,nz,nt) with n^., Uy 
and Hz ranging from to 3 and fit ranging from to 6 for DBW2 f] = 1.22. For DBW2 
f3 = 1.04 we used rix, Uy in the range —2 to 2, in the range to 2 and nt in the range to 
4. These choices gave us a sufficient number of distinct momenta for our non-perturbative 
renormalization. 

Results of the first diagonal element of A^A~^, Z~^Zvv+AAyv+AA^ versus pfj^^t shown 
in Fig. ^Ifor each mj, where the left/right panel is from DBW2 /? = 1.22/1.04 for each mf 
(open symbols) and chiral limit mj = — mres (filled circles). 
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The momentum and mass dependence of the data is expected to be classified into two 
regions as the operator product expansion would suggest. For low momenta, there should 
be significant contributions from hadronic effects. At larger momenta, this effect should be 
suppressed in an approximately power-law manner. For large enough momenta, there is a 
region in which the hadronic effects are negligible and the momentum dependence (running) 
is well described by perturbation theory. As we are working on the lattice, if we work at 
too large a momentum there will be significant contributions due to discretization errors. 
Therefore, the success of the NPR technique requires there to be a window of momenta for 
which contributions from both hadronic effects and discretization errors are small. Our data 
suggest this is the case. Namely, while at very low momenta there is significant momentum 
dependence, in momenta range (~ 2GeV) where we might expect that perturbation theory 
is valid there is only a mild momentum dependence. We use this range in the following to 
compare with the perturbative prediction. 

Figure 1201 shows the results of the off-diagonal element Z~'^ Zvv+aa,vv-aa with the same 
organization as Fig. El For both the case of DBW2 (3 = 1.22 and 1.04, we obtain values 
in the chiral limit that are consistent with zero as a results of linear extrapolation with a 
reliable quality of fit. Figure EH shows examples of mass dependences of Z^"^ Zyy+AAyv-AA 
for some fixed momenta both for f3 = 1.22 (left) and 1.04 (right). Data in this figure show 
linear behavior in confiict with the prediction of 0((mj + mi.cs)^) dependence, implied by the 
discussion in Appendix El This linear behavior should likely be interpreted as an error due 
to spontaneous chiral symmetry breaking. A nonzero value of (qq) could cause systematic 
errors of O ((Aqcd/p)^) and O ((Aqcd/p)'^ ■ ("^/ + ^res)/p), which, for not large enough p, 
contribute as a sizable intercept and a linear term of the quark mass, respectively. This 
interpretation explains that, in Fig. degree of the slopes in both panels and that of the 
intercept in the right panel decrease for larger values of piatt- Another possible source of 
the linear term is the dimension 7 operators which contain one derivative in a four-quark 
operator. 

Other off-diagonal elements Z-'^Zvv+aa,ss-pp, Z-^Zvv+aa,ss+pp and Zg^Zvv+AA,TT 
do not allow a naive chiral extrapolation. For these factors, we plot rrif dependences for 
several points of Pi^^^ in Fig. (221 As can be seen, these elements increase in magnitude as 
quark mass is decreased. While it may seem counterintuitive that these measures of chiral 
symmetry breaking increase for small masses, this effect is well understood to be a result 
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of the poor choice of momenta configurations for the renormahzation condition: to be able 
to ignore hadronic contributions at large momenta, it is necessary that all the momenta in 
the problem are large compared to the hadronic scale, not just the external momenta.^ The 
choice of momenta that we have been forced to use to remain consistent with the perturbative 
matching calculation transfers no momentum through the operator. The contribution of any 
particle that this operator couples to is therefore only suppressed by that particle's mass. 
Practically this is only a problem when the operator couples to a pseudo-Goldstone boson, 
the mass of which goes to zero in the chiral limit. While the presence of these "pion poles" 
greatly complicates any attempt to accurately extract the wrong chirality mixings, in this 
work we will content ourselves with placing a bound on the size of such contributions. As 
such we will extract the mixing coefficients using our heaviest values of the mass and largest 
values of the momentum. In this way we reduce the pion-pole contamination, while, at 
the same time, maximizing the 0{m'j) contributions. As will be demonstrated in the next 
section, the resulting mixing coefficients are small enough that the wrong chirality operators 
can be safely neglected in our calculation of Bk, namely we calculate its renormahzation 
factor as z^^^^'~'^ = Zw+aa/ Z\. 

In Appendix O we summarize the perturbative formulae for the renormahzation group 



running and scheme matching. The authors of Refs. |48l. I49l| calculated a factor absorbing the 
scale and scheme dependence of the renormahzation factor, ^^V^^^^ NLO perturbation 
theory. Using this factor, we convert our results to the renormahzation group independent 
(RGI) value, 

(Nf) = w^//mom(^/. Piatt/ a) ■ ^BK^°^(piatt) • (34) 

This is related to the renormahzation factor at a certain energy scale, /x, in the MS-scheme 
using naive dimensional regularization by: 

^MS NDR(^^^ ^) ^ f^)ZsANf) . (35) 

The functions wri/mom and wjj^ are defined in Appendix [01 and piatt is the magnitude of 
the lattice momentum used in the Green's function defining 2'^^^°^. 



^ This is conventionally phrased as requiring that the external momenta be non-exceptional, i.e. the sum 
of each subset of the external momenta (defined as incoming) must be large. 
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It should be noted that both Eq. EUand Eq.Eil depend upon the number of active flavors. 
While the final result we are aiming at is an estimate of renormalized in full QCD (3 
active flavors), our lattice calculations of the bare value of Bk and the renormalization 
factors in the RI renormalization scheme are performed in the quenched approximation. 
While using Eq. OUand Eq. IHSlwith either Nf = and Nf = 3 - in any combination - would 
seem to be equally valid procedures (just different definitions of the quenched estimate for 
Bx ), in this work we will consistently use Nf = 0. One advantage of this approach is that 
it allows us to compare the observed scale dependence of the Rl-scheme renormalization 
factors versus the perturbative prediction. In this way we are able to study the size of the 
associated systematic error. 

In taking this approach, we must also employ a value of as in the quenched approxima- 
tion. We obtain this using the two-loop formula given in Eq. lC21 with Nf = and = 238 

■' Mo 

MeV. This latter value is gained by taking the value of '"o^^ given by j50| and converting 



it to physical units using ro = 0.5fm, the results of and our quoted value of the lattice 
spacing, as extracted from the chiral limit of the rho meson mass. This is the same approach 
employed in |2l|, where more details can be found. While the quenched coupling constant 



5l| is another possible choice for as, 



obtained from the value of the plaquette and rectangle 
this choice changes the result of by less than 0.1%. 

The comparison between zf^^°^ and Zbj, is shown in Fig. |2Sl for DBW2 p = 1.22 
(left panel) and DBW2 /5 = 1.04 (right panel). One observes that Zb^c is almost scale 
independent for piatt ~ 1- Assuming the perturbation theory at NLO is good enough, the 
remaining small slope of is caused by discretization errors (0((pa)^) effects). To remove 
these errors, we carry out a linear extrapolation in pf^^^ for piatt > 1 and quote the values of 
intercept as Zej^- 

Our final results are 

' 0.9901(36) DBW2 p = 1.22 

(36) 

0.9427(54) DBW2 (3 = 1.04 
Note that it is possible to interpret the small but noticeable slope of the linear fit in the 
right panel of Fig. |2S1 as an error of the perturbation theory at lower momentum of ~ 2 
GeV. Taking the constant fit instead, we obtain ~ 1% larger value for (3 = 1.04. We may 
compare numbers in Eq. EHl against one- loop lattice perturbation theory calculation which 



^MSNDR^^^^Q^ ^ = 2 GeV) 
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has been done in Ref. jsij. Using measured values of the plaquette and rectangle as input, 
the perturbative renormalization factors are 



^MSNDR(^^_Q^ ^ = 2GeV) 



0.9775 DBW2 /3 = 1.22 



(37) 



0.9493 DBW2 /3 = 1.04 

both of which lie close to the non-perturbative values. We find it reassuring that these two 
quite different methods lead to results agreeing to better than 2% 



VI. KAON S-PARAMETER 



A. Chiral behavior of K—K matrix element 



Before presenting our results for Bx-, it is important to check the chiral behavior of the 
K-K matrix element of Q^^^^'^\ We calculate the three point correlation function with 
degenerate quarks and find a suitable plateau for ti <^ t <^ t2 to extract the desired matrix 
element: 
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X 2mps. 



(38) 



ti<i<i2 



As mentioned earlier, the sink and source locations were chosen to be (^1,^2) = (7,41) for 
DBW2 (3 = 1.22 and (5,27) for (3 = 1.04 and Wilson /? = 6.0. Results from a constant fit 
to the plateau of the matrix element for each m/ in the fitting ranges 19 < t < 29 (DBW2 
p = 1.22) and 14 < t < 17 (DBW2 p = 1.04 and Wilson p = 6.0) are hsted in Table IXTTl 
If a single intermediate state contributes to the matrix element in the numerator of the 
right-hand-side of Eq. IHHl this quantity will be independent of the intermediate time t. This 
is demonstrated in Figs. El and 123 which show this quantity as a function of t for each of 
the masses analyzed. Both graphs show an apparent plateau region which could be as large 
at 17 time units for /? = 1.22 and 10 time units for /? = 1.04. It seems likely that in the 
fitting range chosen contamination for excited states should be below a few percent. 

In quenched chiral perturbation theory this matrix element, expanded in powers of mpg 
up to O(mpg), is given by 
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where, following the discussion in Section IIIH we neglect both the quenched chiral log and 
the chiral log terms in mp^. Since the expansion starts at 0(m|g), this matrix element 
vanishes in the chiral limit, mpg = 0. Another characteristic of Eq. [HHlis that the ratio of 
the coefficients of the leading term and the chiral-log term is determined solely by /, the 
decay constant in the chiral limit. It is interesting to examine our data in light of these 
expectations from quenched chiral perturbation theory. For this purpose, we carried out the 
two-parameter fit to Eq. EHl using for / the product of chiral limit value /o and listed 

"2'' PS^ in lattice units and the fitting curve from 
Eq. ISnifor DBW2 /3 = 1.22 (left panel) and /? = 1.04 (right panel). We also used the fitting 
functions 



in Table H Figure Ei shows (PS Q^^^ 



/PS 



Q 



(AS=2) 



PS 



1 



6 



m 



PS 



(40) 



and 
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to examine possible explicit chiral symmetry breaking effects through the magnitude of 
ao and to compare the result for as/ai to the chiral perturbation theory prediction of 
-6/(47r/)2. As listed in Table IXTTTl all three fitting functions in Eqs. El EOl and ED with 
A^ = 1 GeV fit our data equally well. Results for a\ and 02 are consistent among the fits, and 
ao in Eq. HUl is consistent with zero. The latter agrees with previous quenched domain-wall 



)(A5=2) 



PS) vanishes in the chiral limit, mpg 



0, 



fermion calculations that showed <^PS |( 

or m/ = — mres, to good accuracy [2], The same is true of our recent calculation of 
Bk in the two-fiavor theory [2^. Furthermore, 03/01 from Eq. 1^ reproduces the analytic 
result — 6/(47r f)2 fairly well, as did our previous calculation using the Wilson gauge action 
at /? = 6.0 l21]. This was not the case in Ref. [s^ though in that study the authors did 



not examine (PS 



=2) 



,(AS=2) 



PS) / (0 \P\ PS)^ where 



PSy directly, but a ratio (^PS Q 
P is the pseudoscalar density. While it is reassuring that our data fits standard quenched 
chiral perturbation theory so well, it should be kept in mind that the analysis presented 
here includes quite heavy pseudo-scalar masses which may lie above the region where chiral 
perturbation theory is valid and has neglected quenched chiral logarithms and finite volume 
effects which may distort the lightest mass points. 
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B. Results for Bk 



Let us now discuss our results for the kaon S-parameter, Bx defined in Eq. ^ Following 
our earlier conventions for the decay constant and pseudoscalar mass, we will use the notation 
Bps for this AS = 2 amplitude evaluated for pseudoscalar states with a general meson mass, 
mps- The parameter B^ will be used when mps = m^: B^ = Bps{mx)- In the lattice 
calculation of Bp^ we deal with the same three-point correlation function as in the previous 
subsection, but in a ratio with two factors of the wall-point correlation function C^^, 
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(42) 



Here we suppress the appropriate (latt) superscript, -Bpg**-* Bps- 

Plateaus for this ratio for each value of are shown in Fig. EZI for DBW2 j3 = 1.22 
and in Fig. |2Hlfor j3 = 1.04 (the quark mass increases from bottom to top in each figure). 
The solid and dashed lines in each plot indicate the fitting range used (19 < t < 29 for 
f3 = 1.22 and 14 < t < 17 for j3 = 1.04) and the results for a constant fit that are also 
listed in Table IXIVI As is suggested by Figs. EH and EH these fitting ranges are chosen 
quite conservatively and could likely be made larger without significant contamination from 
higher mass, excited states. In fact, choosing a larger fitting range yields consistent results. 
For example, for our lightest masses, increasing the fitting range to 17 < t < 31 for f3 = 1.22 
decreased the result by 1% while for (3 = 1.04 the enlarging the fitting range to 12 < t < 19 
increased the result by 2.5%, both within one standard deviation of the results quoted in 
Table IXlV 



In Ref. j2l|, we chose a different method to determine -Bps, computing separately the 
numerator and denominator of Eq. Upland then evaluating their ratio. In the present case. 



that method and the one used here give results which agree within statistical errors. Another 
variant of our method replaces the quantity /pg formally contained in the denominator of 
Eq. E21 with the alternative /pg* which is computed in Section HVl However, as can be seen in 
Fig. the difference between the two is always less than two percent, even for the lightest 
ruf, and usually smaller than one percent. 
The counterpart of Eqs. EHl and El] for Bk is 
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and 



2 

Bps = ^0 + ^imls + 6"^ps In (44) 

X 



respectively |4^ . For degenerate quarks, these have the same form as in the theory with sea 
quarks js^, i.e. there are no quenched chiral logarithms in Bps because of the cancellation of 
such terms between numerator and denominator in Eq. 021 Values of the fitting parameters 
for these functions are given in Table IXVl As seen in this table, we find both fits are 
equivalent, with the /3 = 1.04 case showing the closest agreement. 

We also constructed Bps from (fS Q^^'^^^^ PS^ calculated as described in the pre- 

(2) 

vious subsection and /pg obtained from a constant fit to the plateau in t for the ratio 
Cp^^(t,t2)Cp^^{t,ti)/C^^(t2,ti). While the central value of the result changes by less than 
0.2%, the jackknife error on the ratio increases by ~ 70% compared to the jackknife error 
coming directly from the use of Eq. 

In Fig.|2ni we plot bare values of -Bps versus mpg (DBW2 (3 = 1.22 and 1.04 for the left and 
right panels, respectively). The solid and dashed curves denote fits to Eqs. 1431 and 1^ which, 
for DBW2 /? = 1.22, somewhat differ in contrast with the case of the (PS Q^^^=^^ PS) 
matrix element where they were almost on top of each other. We found that this difference 
does not depend on the choice of fitting range. Repeating the same analysis for both f3 
values using the larger fitting ranges described above did not change the situation. Since 
we extract Bk from an interpolation of Bps to the kaon mass (around the data point for 
second lightest mass), the choice of fitting function makes little difference in the value of 
Bk, as discussed below. In the absence of a compelling reason to choose one over the other, 
we pick Eq. 031 to be consistent with chiral perturbation theory. 

Interpolation to the physical point, mps = mx = 495 MeV, yields the lattice values of 
^(latt) jjg^g^ jj-^ ^Yie last column of Table IXVI and indicated by the open symbols in Fig. (211 
Though we use Eq. 031 to obtain B^^^^^ in the rest of this article, the difference from using 
Eq. Elis always less than 1%. After multiplying Bps by ^^'^{Nf = 0, /i = 2 GeV) in 
Eq. |2ni we can directly compare the renormalized values from each ensemble as shown in 
Fig. EDI where filled symbols denote DBW2 /? = 1.22 (circles) and P = 1.04 (squares) and 
open diamonds, Wilson (3 = 6.0. Fitted curves corresponding to Eq. E31for DBW2 /3 = 1.22 
(solid) and (3 = 1.04 (dashed) are also shown in the figure. Since the points do not lie along 
identical curves, there are evidently lattice spacing errors remaining in our determination of 
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Bk. 

As was discussed in Section a possible systematic error that may contaminate is 
mixing with wrong chirahty operators through renormahzation: 



bt' = bk+y: ^VV+AA, ./Zl-Br\ (45) 

i 

where the S-parameters for the wrong chirahty operators, Oi = Ovv-aa, Oss±pp, and 
are defined as 



p(iatt) _ \PS|C»|PS^ 

- 8.2 ^2 • (46) 

In Section El we pointed out the large, 0(m^gj suppression of this contamination. Here, we 
press our point by a numerical demonstration. 

For DBW2 (3 = 1.04, we have calculated all of these 5-parameters^ following the same 
methods used for -Bps- The results are listed in Table IX VII and shown in Fig. |^ The 
magnitudes of the i?-parameters for these wrong chirahty operators are less than two or- 
ders of magnitude larger than -Bps, even for quark masses of ~ This makes their 
effects at ms/2 very small, given the further O(m^gg) suppression present in the mixing co- 
efficients. The difficulties in accurately determining these miniscule mixings were outlined 
in Section namely we can not measure them accurately with our current techniques. As 
a gross overestimate of the size of these mixings, we measured their values at the largest 
momentum p^^^^ = 2.4674 and the heaviest quark mass fnj = 0.05 and find contamination 
from each Oi is no larger than 0.01. Moreover, cancellation between the wrong chirahty 
terms in Eq. ^1 likely makes the net contamination even smaller. Thus, we conclude that 
in our determination of Bk, the contributions from the wrong chirahty operators are well 
below our statistical error. 

Results for B^ ^^^(/i = 2 GeV) = (iV/ = 0, ^ = 2 GeV) ■ -BJ''"^ and Bk = Zb^ {Nf = 
0) ■ 5^^**'' are collected in Table IX VIII where we enumerate perturbatively (PR) as well 
as non-perturbatively (NPR) renormalized values. Results of a linear extrapolation and a 
constant fit to the continuum limit for each quantity are listed in the first two rows. As 
mentioned in Section El while B^^ is almost independent of our choice for A^, Bk is 



^ Due to their relevance for beyond-the-standard-model physics, these results may also be useful for future 
studies. 
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significantly sensitive. For that reason, we focus on the result for the former quantity in the 
following. 

Our results for ^°^(/i = 2 GeV) are shown in Fig. |221as a function of aL 
results from Ref. [2l| and the results obtained by the CP-PACS collaboration 
domain-wall fermions with parameters similar to ours (a~^ = 2.87 GeV, 24^ x 60 sites, 
Lg = 16 and = 1.88 GeV, 16^ x 40 sites, Lg = 16). The main difference from our 
calculation is their use of the Iwasaki gauge action and perturbative renormalization of 
i?pg**\ At a^^ ~ 2 GeV, results from the two collaborations differ by roughly two standard 
deviations: ndr^^ _ 2 GeV) = 0.564(14) (CP-PACS) and 0.532(11) (Wilson f3 = 6.0), 
0.524(11) (DBW2 j3 = 1.04). The two results are even more consistent at the smaller lattice 
spacing. 

The discussion of the continuum extrapolation of the J-parameter and the decay constants 
in previous sections is valid here as well. While the result of constant fit in Table IX VIII 
is quite acceptable with a x^/dof of 1.8, we use the linear extrapolation, shown in Fig. 
to obtain our final result because we have no a priori reason to expect the term to be 
absent. This linear extrapolation is done by connecting our two data points and the error 
for the continuum limit is obtained by quadrature. We obtain the final result: 

B^ ^°^(/i = 2 GeV) = 0.563(21)(39), (47) 

where first error is statistical and the second systematic, which is discussed in the next 
subsection. 

Finally we give a result for -Bps, evaluated in the chiral limit and then extrapolated to the 
continuum limit: i?^^^°^(mps = 0, fi = 2 GeV) = 0.289(18), where we do not attempt to 
determine the systematic errors for this value because of the sizable uncertainties associated 
with evaluating the chiral limit from our data. As is displayed in Fig. 123 we again see 
relatively mild dependence on the lattice spacing. This result is based on the chiral fits using 
the known chiral logarithm given in Eq. 1321 and tabulated in Table IXVI Here we do not 
compare with the CP-PACS result for this quantity because their use of an un-constrained 
fit to the chiral logarithm, which we are able to avoid, introduces large uncertainties in the 
chiral limit. 
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C. Estimate of systematic errors 

All of the statistical errors obtained earlier in this paper have been assigned using standard 
jackknife procedures and should reflect the variations that would be seen if our Monte Carlo 
calculations were simply repeated and analyzed in an identical fashion. (Of course, one 
should recognize that these errors, obtained from a finite data set are themselves subject 
to error.) Less certain but equally important are the errors in our results that come from 
systematic limitations in the calculation. These may be crudely divided into two types. 
First and easiest to determine are those associated with our methods of analysis. Our 
choice of plateau region, the procedure for chiral extrapolation or interpolation and our 
method for taking the continuum limit are good examples. Here the calculation should 
show consistency with the theoretical ideas being used in the analysis and the variation in 
the result between different approaches should indicate the level of systematic error. Of 
course, if the theoretical framework describes the data poorly or contains many parameters, 
a reliable estimate of systematic errors may not be possible. 

The second type of error reflects ingredients which are missing in the calculation. If 
only a single volume or lattice spacing is used, the errors associated with finite volume 
or finite lattice spacing cannot be determined from the calculation at hand. Similarly the 
errors induced by the quenched approximation cannot be known if no full-QCD calculations 
have been performed. While one may "estimate" an expected error by comparing with 
more extensive calculations of other quantities, such estimates often reduce to an exercise 
in wishful thinking. However, for the case of Bk there are now many results reported from 
other Bk calculations which either independently, or by comparison with the work presented 
here, provide reasonably direct information about all of the important sources of error. 

In the discussion to follow and the final results quoted we attempt to estimate the size 
of these systematic effects based on the calculations presented in this paper and the results 
of other work. These "systematic errors" are not intended to be upper bounds on the size 
of these systematic effects but an estimate of their likely size. Thus, in performing such 
estimates we will not add cautionary inflation factors as would be appropriate if we were 
attempting to deduce reliable upper bounds on these errors. Rather we believe that we 
can extract the greatest value from this calculation by attempting to directly determine 
the suggested size of these effects. Attempting to determine the errors on these estimates 
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(finding "errors on errors") or to establish reliable upper bounds on the size of these effects 
is beyond the scope of the present work. Each possibly important source of systematic error 
will now be discussed in turn and the results summarized in Table IXVIIII 

Extraction of K — K matrix element. As was discussed in Section IVl Bl there is system- 
atic uncertainty inherent in our method of determining the ratio Bps defined by Eq. |321 
associated with the choice of fitting region. Specifically, if the fitting range includes times 
too close to the source or sink, the resulting value for Bps may receive contamination from 
excited states. This was studied quantitatively by comparing two choices of fitting region 
where variation in the result for -Bps, on the order of the statistical uncertainty was seen. In 
this situation, our choice of fitting range determines the character of this error. If we had 
chosen a large fitting range, risking such excited state contamination, we would see a smaller 
statistical error (reflecting the larger number of points in our fit) but a larger systematic 
error. The systematic error would be determined by a comparison with the smaller fitting 
range and would likely be dominated by the statistical error from the smaller fitting range. 
In the approach we have taken, using the safer, smaller fitting range, the error is essentially 
statistical since we are well away from a region where excited state contamination might be 
expected. As might be deduced from Fig. |211 the data shows so little time dependence, that 
it is not possible to reliably extract the mass of a possible excited state from the relevant 
correlation functions. Instead, we estimate this possible contamination by evaluating e"^™'* 
for t corresponding to the 11 lattice spacings, the distance between our measurements and 
the source. If we choose Am = 1 GeV as the gap between our K meson and the first ex- 
cited state with the same quantum numbers, this suggests an upper bound on this possible 
percentage contamination of 3%. 

Determination of fx- Since B^ is the ratio of a. K — K matrix element divided by f]^, the 
systematic error in determining the kaon decay constant must enter our result for B^- This 
was discussed in Section HVl where, by carefully comparing statistically correlated quantities, 
we were able to recognize a systematic difference between different methods for determining 
fx, possibly caused by unsuppressed near zero modes. These differences were on the level 
of 1% for second to the lightest mass which has the greatest effect on the determination of 
-Bps at the Kaon mass. Thus, this source of error is listed in Table IXVlIII as a 1% effect. 

Kaon mass. Although Bk is dimensionless, it is obtained by interpolation to the point 
mps = rriK and therefore depends on our choice for the K meson mass. While we have 
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determined the Kaon mass directly in lattice units quite accurately (~ 1%, see Tables HVl and 
\V\ there is further uncertainty in determining the lattice scale in physical units, especially in 
a quenched calculation. Following past practice, we have determined the lattice scales given 
in Table IVIII from the p mass. However, the fact that the p meson is a stable state in this 
calculation but is an unstable particle in Nature with a width to mass ratio of 20%, suggests 
that this may introduce significant systematic errors. The dimensionful decay constants 
and fx provide alternative values for the lattice spacing. The discrepancies between 
their continuum limits (using rUp to set the lattice scale) given in Table IXII and experiment 
provides a simple ~ 6% estimate of this source of systematic error in the choice of value 
for rriK- Referring to the dependence of Bk on the input kaon mass shown in Eq. 1^ or 
extracted from Tables HV ] IVl and IXIV| we conclude that a 6% error in tjik propagates into 
a 3% error in B^- 

We can also use the static quark potential to set the lattice scale. This was done for the 
quenched lattice configurations studied here by Hashimoto and Izubuchi in Ref. 

Q- The 

comparison between the lattice scale determined from rrip and that implied by a choice for 
the Sommer scale of tq = 0.5 GeV is given in their Table 1, showing agreement on the 6-9% 
level, roughly consistent with our 6% estimate. As mentioned above, much of the difficulty in 
determining the lattice scale from experiment comes from the quenched approximation and 
hence may already be represented in the error associated with the quenched approximation 
discussed below. However, we have adopted the conservative approach of listing the effects 
on Bk of this resulting uncertainly in rriK as a separate error. 

Operator normalization (NPR). The non-perturbative renormalization of the left-left op- 
erator Q^^^^"^^ as it appears in Eq. ^provides a case where we expect the RI/MOM pro- 
cedure of the Rome-Southampton group to be particularly accurate. The somewhat less 
precise wave function renormalization constant Zq cancels in this ratio and, as can be seen 
from Figs. |2S1 there is a large kinematic region, 1 < piatt? free of infrared QCD effects. 
We believe that the principle source of systematic error in the determination of Zbj^ is the 
presence of lattice spacing errors in this region. As discussed in Section El we attempt to 
remove these errors by identifying and subtracting an (piattc^)^ term. This introduces a 1% 
change which we will adopt as an estimate of the systematic error in this determination of 

Zbk- 

Operator normalization (PT). As is discussed in Section |3 a final perturbative step is 
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needed to convert Bk defined in the RI/MOM scheme to the more conventional, pertur- 
batively defined MS NDR scheme. This correction factor is 1 at tree level while the NLO 
correction, 0{as), introduces a 1.3% change. In order to provide a proper estimate of the 
omitted 0(a|) term we would need either a two loop result, which is not available, or a 
step-scaling connection between operators normalized at the present lattice scale and those 
normalized at a much finer scale where one-loop perturbation theory will be more accu- 
rate ^53^. Lacking both of these alternatives, we will make the conservative estimate that 
this continuum change of scheme factor at two-loops is literally q;| = 0.04 or a 4% correction, 
three times as large as the admittedly small one-loop correction. Note, there is a further 
uncertainty in our one-loop correction coming from our choice for since we must again 
determine a zero-fiavor quantity from experiment. We will not add a further systematic un- 
certainty from this source, relying instead on the error associated with quenching discussed 
below to include this effect. 

Mixing with wrong chirality operators. This was discussed at length in Section IVI Bl Our 
attempt to estimate the chirality violating mixing coefficients numerically give a potential 
error of order 0.01. However, since this numerical result represents an upper bound on 
quantities too small to be computed with our present resources and there are good theoretical 
arguments that the mixing coefficients should be on the order of 10~^, a more accurate upper 
estimate for this sort of error is 10~^ or 0.02%. 

Finite volume. The calculations described in this paper were performed using a single 
relatively small physical volume of approximately 1.6 Fermi on a side. Thus, we cannot 
estimate the errors associated with this choice of volume from the results presented here. 
Instead, we use the results of Ref. bc| which provides values for both 1.7 and 2.6 Fermi 
volumes, seeing a ^ 2% increase in the result for Bk from the larger volume. We will 
interpret this difference as the finite volume error in the result presented here. 

Degenerate quarks in the Kaon. In contrast to Nature, the K meson state studied in this 
paper is composed of two degenerate quarks each with mass approximately one-half that of 
the strange quark. It is known empirically that to a good approximation many quantities, 
Bx included, depend only on the sum of the constituent quark masses. In a recent RBC 
collaboration paper |3l| the effect of non-degenerate quarks was resolved explicitly and found 
to cause a downward shift of B^ of approximately 3%. Since that was a full QCD calculation 
with two fiavors of dynamical quarks at a single lattice spacing it can not be precisely related 
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to the calculation presented here and we include this shift as a further systematic ±3% error. 
(Note, the presence of quenched chiral logarithms prevents this question from being studied 
in the quenched approximation.) 

Continuum extrapolation. As can be seen from Fig. |22l our values for at a^^ ~ 2 
and 3 GeV, show relatively mild lattice spacing dependence. Since domain wall fermions 
are expected to show finite lattice spacing errors of order a^, the linear extrapolation shown 
in that figure should provide a good estimate of the continuum limit. The largest possible 
source of systematic error in such an extrapolation that involves only two points is the 
additional term of higher order lattice spacing. Of course, with results for only two values of 
lattice scale we cannot determine such term. To get a rough idea of its size, we will assume 
that the O(a^) (known) and O(a^) (unknown) terms are the first and second terms in a 
geometric series in a^. (Note the absence of 0{a^) errors results from the chiral symmetry 
of domain wall fermions.) This approach implies a higher order correction of less than 1% 
and changes the extrapolated value of Bk by 0.0011 or 0.2%. Note, this is an estimate of 
the systematic error in our result for Bk coming from the continuum extrapolation. There 
is also a contribution to the overall statistical error coming from this extrapolation which 
has been incorporated using standard error propagation. 

Omission of quark loops. The most important source of systematic error in the results 
reported here comes from our use of the quenched approximation. While the present 
calculation is entirely quenched we can compare this quenched result with a similar do- 
main wall calculation recently performed on a 16^ x 32 volume with lattice spacing of 
~ 1.7 GeV and two flavors of dynamical fermions This calculation gave a 

value for i?^^(2GeV) = 0.495(18). Because this result was obtained on a coarser lattice 
(a~^ = 1.691(53) GeV) than the present one, we should extrapolate our present result to this 
larger lattice spacing. Applying the same linear in form used to determine the continuum 
limit, the present quenched calculation predicts the value i?^^(2GeV) = 0.512(19). The 
comparison can be seen graphically in Fig. 1221 where both the results of the present work 
and the earlier Nf = 2 result are shown. 

Since the statistical errors on these two results and their difference are all the same, 
~ 0.020 there is no evidence for a systematic shift caused by including the effects of quark 
loops. Because the fuU-QCD calculation we are using for comparison contains only two light 
quarks, not the three quarks present in Nature, we attempt to account for this discrepancy 
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by increasing this error estimate by a factor of 3/2 and include a symmetrical error of ±0.030 
in Table IXVIIII to represent the systematic error coming from our omission of quark loops 
in the present calculation. 

The implication and interpretation of Table IXVIIII are discussed in the conclusion. 

VII. CONCLUSIONS 

We have presented a study of quantities related to the physics of the light quarks u, d 
and s using, primarily, results from the DBW2 gauge action with gauge couplings, P = 
1.22 and 1.04 (a^^ ^ 2 and 3 GeV, respectively). Due to the combination of domain- 
wall fermions and DBW2 gauge action, the residual quark mass mres is small compared to 
the input quark masses thus ensuring that important chiral symmetry properties in these 
simulations are intact. Infrequent tunneling between topological sectors as a ^ is made 
worse by the DBW2 gauge action since small dislocations that aid the tunneling process are 
suppressed These are precisely the configurations that lead to relatively large chiral 

symmetry breaking effects. To avoid the problem of an excessively large number of sweeps 
between pseudo-independent configurations, for the finer lattice [jS = 1.22) we adopted 
a new strategy of using many different initial configurations generated with the Wilson 
gauge action which allows more frequent tunneling. Thus, the entire ensemble of DBW2 
configurations reflected the initial, more physical, distribution of topological charge of the 
Wilson lattices. It is worth emphasizing that this is an algorithmic problem, not a deficiency 
special to the DBW2 action: as a — > all actions updated with a small step algorithm will 
tunnel less and less frequently. 

The present calculation contains considerable information about the size of possible sys- 
tematic errors. Earlier CP-PACS results on multiple volumes and RBC Nf = 2 QCD 
results j^^, including two flavors of dynamical quark loops, provide further information 
about possible errors associated with finite volume and quenching. We will now combine 
this information, summarized in Table IXVlIII to deduce a final value for B^- We begin with 
the direct result of this calculation, given in Table IXVTTl B^^^^{fi = 2 GeV) = 0.563(21), 
where the quoted error is statistical. We incorporate the discussion of systematic errors in 
the previous section as followings. First we combine in quadrature all systematic errors, 
with the exception of quenching, into a single systematic error of ±0.039. For clarity, we 
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then quote the quenching error separately as ±0.030. This estimate does not come from a 
direct comparison of a quenched and full QCD calculation with the proper number of flavors 
performed at the same lattice spacing. Rather, as discussed above, it is a bound on the 
possible error coming from a comparison with a 2-fiavor calculation performed on a coarser 
lattice and the quenched calculation reported here extrapolated to the same coarser lattice 
spacing of a^^ ~ 1.7 GeV. Since the quenched and 2-fiavor numbers agreed within errors, 
we are unable to see an effect from the inclusion of 2-flavors of dynamical quarks and use 
the difference of these two numbers as an estimate of the quenching error in the present 
calculation. Thus, our final result is 

^MSNDR^^ = 2 GeV) = 0.563(21)(39)(30), (48) 

where the first error is statistical, the second systematic (excluding quenching effects) and 
the third represents the quenching error. 

While we believe that it is important to attempt to assess the possible systematic errors 
in our calculation of Bk-, this exercise should be viewed with considerable skepticism. In 
addition to the obvious limitations in many of these estimates, our treatment also ignores 
possible couplings between the different types of errors. Certainly the errors coming from 
finite volume and finite lattice spacing are based on quenched calculations and could be 
significantly different for a full QCD calculation with light dynamical quarks. Were the 
finite volume correction for Bx to double for the case of dynamical quarks that 6% effect 
alone would equal our estimate of the systematic error. 

There are similar concerns regarding the finite lattice spacing and quenching errors. One 
might argue that the O(a^) errors present in the dynamical calculation of Ref. ^jj which we 
use to estimate our quenching error are determined by coefficients of dimension-6 operators 
in the Symanzik effective Lagrangian that describes this lattice theory in the continuum. In 
the treatment presented here, we are assuming that these coefficients can be evaluated with 
reasonable accuracy in the quenched approximation. However, the short distance effects 
which are summarized by these coefficients are quite different for full and quenched QCD 
because the short-distance screening effects of the quark loops in full QCD are potentially 
large. Both of these topics require more study and understanding before the above estimates 
of systematic errors can be treated as trustworthy. 

Our results for Bx are smaller than those reported previously from calculations with Wil- 



son or Clover fermions 
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56 



57 



39 



58| and staggered fermions jlO'] by roughly ten 



percent. Given the uncertainties associated with chiral symmetry in the Wilson or Clover 
fermion calculation, it is most appropriate to compare with the quenched staggered calcu- 
lation of Ref. Q which obtains = 2 GeV) = 0.628(42), roughly two standard 
deviations larger than our value quoted in Eq. HHl While this discrepancy could be caused 
by an unlikely statistical fluctuation, we speculate that it is more likely the result of system- 
atic effects, perhaps the difficulty of performing the staggered continuum extrapolation in 
the presence of large scaling violations. We are encouraged that this result and the previous 
quenched domain-wall fermion calculation of CP-PACS with similar physical volumes but 
different gauge actions agree in the continuum limit as indicated in Fig. |221 Likewise, there 
is nice agreement between our ^ 2 GeV calculation performed using the Wilson and 
DBW2 gauge actions. 

Note, that Eq. HHl presents a hybrid result for in which a numerical estimate of 
quenching errors performed at a single lattice spacing is combined with a continuum limit 
taken from a quenched calculation. Such a result depends on the assumption that the effects 
of quenching (here described as a 4% uncertainty) result in a correction of similar scale on 
the effect of taking the quenched continuum limit (which causes a 7% increase). That is the 
quenching error in the shift caused by taking the continuum limit might be estimated as 4% 
of 7% or 0.3 %. Under this hypothesis, the value quoted in Eq. HSl adequately accounts for 
both the effects of quenching and the continuum limit. Of course, a full QCD calculation at 
a variety of lattice spacings is ultimately needed to remove the uncertainly associated this 
assumption, a task now being undertaken with the next generation of computers. 
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APPENDIX A: OPERATOR MIXING 

In this section we describe the size of the mixing, under renormahzation, between op- 
erators in different chiral multiplets. A basis for the complete set of operators in question 
is Ovv±AA, Oss±pp and Ott defined in Eqs. OilOHl If chiral symmetry was only (softly) 
broken by the fermion mass, then these operators form three distinct sets, which do not mix 
under renormahzation. In particular, the operator of interest here, Ow+aa would renormal- 
ize multiplicatively. However, when using the domain-wall fermion formalism, each operator 
mixes with all the others under renormahzation due to the (small) explicit breaking of chiral 
symmetry. As the relevant matrix elements of the wrong-chirality operators are much larger 
than the one we are interested in, it is important to have a method of estimating these 
mixings. 

Under certain reasonable assumptions, the size of the mixing coefficients can be estimated 
in terms of the residual mass. A framework for understanding how the explicit breaking of 
chiral symmetry is manifest in the associated QCD Lagrangian is outlined in Ref. (2fll |. 
Using the notation and conventions of Refs. j^,!^' introduce an additional term in the 
domain-wall fermion action which reads 

Sn = -Y. ^'^,L,/2-iPl (^^^ - l) ^x.,L./2 + '^x,LJ2Pr - 1) ^x,L./2-l] • (Al) 

X 

The modified action possesses an exact, spurionic, symmetry under SU^Nf)^ ® SU{Nf)ji 
transformations if the constant unitary, Nf x Nf matrix Q transforms as 

n-^Unnul, (A2) 

while the standard domain- wall fermion action is recovered in the limit ^ 1. Assuming 
that the effects of explicit chiral symmetry breaking are local, we can then analyze the form 
these effects may take in a low energy effective Lagrangian describing the physics of the 
modes bound to the two walls by studying the operators allowed by the generalized chiral 
symmetry of Eq IA21 An instructive example is the leading order modifications to the QCD 
Lagrangian itself: at the lowest order in the lattice spacing, the most general form of the 
relevant term that can be added to the QCD Lagrangian density modifies the mass term to 
read 

Zmmfqq + c qO) PrQ + qVLPiq , (A3) 
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where c is a constant with mass dimension one. Setting = 1 it is easy to recognize that 
c = Zmmj-cs- Since m^es contains a single factor of Q, we should expect its suppression to 
correspond to a single propagation from one wall to the other in the fifth dimension: the 
minimum propagation needed to encounter the factor of Q, introduced at the mid-point 
s ^ Ls/2 in Eq. El 

For smooth gauge fields, for which the domain-wall fermion mechanism is working well, 



we can simply associate a suppression factor of O(mrcs) for each factor of Q 20|,|27[. Noting, 
from Eq. lA21 that the effect of Q (Q^) is to "flip" the chirality of a fermion into the opposite, 
we may motivate this result in a physically intuitive way: for a fermion to flip chirality due 
to the explicit chiral symmetry breaking of domain-wall fermions, the two walls must be 
connected through the bulk of the fifth dimension. Each such trip through the bulk comes 
with a suppression factor due to the small overlap between the wavefunctions of the quarks 
bound to the two walls. The size of this suppression factor can be estimated by simply 
measuring mj-es, which, as arg ued above, is associated with a single factor of Q. 

As pointed out in Ref . , this analysis can be more complicated away from the smooth 
gauge field limit. In this case the Hermitian Wilson Dirac operator can have a significant 
number of near-zero modes in the range of (negative) mass corresponding to the domain- 
wall height used in simulations. This operator is closely related to a transfer matrix that 
may be constructed describing propagation in the fifth dimension ^| and such zero-modes 
correspond to eigenmodes of this transfer matrix with near unit eigenvalues, implying un- 
suppressed propagation in the fifth dimension. Such unsuppressed propagation in the fifth 
dimension may invalidate the argument which associates a simple factor of O(mres) with 
each factor of Q. 

However, this transfer matrix description provides a more refined language that can be 
used to make a similar analysis. We expect that the residual mass mres arises from modes 
of the transfer matrix that are of two types: relatively plentiful modes (extended in four- 
dimensions) with transfer matrix eigenvalue substantially below one, < e^'^'^, and the rare 
modes (localized in four dimensions) mentioned above with transfer matrix eigenvalues close 
to unity. (Here Ac is the "mobility edge" of Golterman and Shamir jo^.) Contributions from 
the former are suppressed exponentially oc e'^"^" while those from the latter are only power 
law suppressed, oc 1/Ls, but are further suppressed because such localized modes are rare. 

As demonstrated below, the operator mixing of interest here involves changes of chirality 
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by two units and requires two powers of the matrix Q. In such a circumstance, these two 
different types of modes may contribute differently. Extended modes with exponentially 
suppressed propagation will naturally contribute to operators changing chirality by two units 



as m'^cs with each transversal of the fifth dimension introducing a factor oc e ^''^^ . However, 



it may be possible that a single localized mode with transfer matrix eigenvalue close to one 
could support these two transversals with only the suppression implied by the presence of 
that one, relatively rare mode. Never-the-less, a more detailed analysis, to be given in a 
later paper, shows that this is inconsistent with either Fermi statistics or baryon number 
conservation. Such chirality-two changing processes require either two distinct modes for 
the propagation of two quarks of the same flavor or for a particle and an anti-particle. 

Given the preceding discussion, we may reformulate the statement that the mixings 
between Ovv+aa and the wrong chirality operators is suppressed by a factor of 0(m^gj, as 
the statement that we must flip two left-handed quarks into right-handed quarks (or vice- 
versa) to move between these two sets of operators. To explicitly derive this statement, we 
re-write the basis operators in terms of the left- and right-handed components of the quark 
field. Using the chiral representation of the gamma matrices. 



Qr 










1 




, 75 = 











-1 



(A4) 



with a = (1, —i<Ji^2,3) and a = (1, i<Ji^2,3)- Up to overall numeric factors, these operators 
now read 



Ovv+AA OC SLO-f^dLSLO-f^dL + S rW f,d rS rW f,d r 

Ovv-AA OC 'sicr^dLSRa^dR 

OsS-PP OC SLdRSRdL 

Oss+pp OC SLdRSidR + SRdLSRdi 

Ott oc -SRA^''dL-SRA^''dL+-SLA^''dR-SLA^''dR, 



(A5) 
(A6) 
(A7) 
(A8) 
(A9) 



where in Eq. IA9I we have introduced the notation 



(AlO) 
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A'"' 
A^''' 

It is now simple to deduce that to move between Ovv+aa and the operators in Eqs. IA6I - 
IA9I requires two flips of chirality. 



APPENDIX B: RENORMALIZATION CONDITION 



The basic building block for the construction of the quantities needed for the renormal- 
ization calculation is the quark propagator from a single point source ( in this case situated 
at the origin of our co-ordinate system) to a point sink in a fixed gauge, S{x\0). This 
propagator is Fourier transformed into momentum space on the sink co-ordinate 

X 

The amputated vertex functions needed can then be easily constructed in terms of 



-1 



a/3 



and 



G'ip) 



75 (Gip) {Gip))-'yj, 



af3 



(B2) 



(B3) 



where angled brackets represent the average over gauge configurations. For example, the 
amputated Green's functions of the flavor non-singlet bilinear operators 



{u{p) uTd d{p)) 



amp ; 



(B4) 



where F represents the particular gamma matrix (either vector, axial-vector, scalar, pseu- 
doscalar and tensor) may be written as (V^), {A^), {S), {V) and {Tf^u), where 



V^.(p)a/3 

A.^{p)a(3 
^{P)al3 
V{p)ap 



[G'{phAp)U 
[g'{pHi,g{p)U 

[G'{p)G{p)U 

[g'{pHG{p)U 
[g'{pKMp)U- 



(B5) 
(B6) 
(B7) 
(B8) 
(B9) 
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Amputated Green's functions for the relevant operators in this paper are 



p(2)latt/ X 



^(4)latt / N 
VV+AA\P)al3;-yS 



r 



(4)latt 
VV-AA 

T(4)latt 



/3;7<5 



SS-PP\P)a 



/3;7<5 



T(4)latt 
55+PP 

p(4)latt 

i rprp 



{P)al3;-y5 



i.P)a 



{A)ap (BIO) 

(V/^,a/3V^,75 — "^fi,Q/3-^/i,75 — V^,a5V/_j,'y/3 + ^^,05-^^1,7/3) (B12) 

{Sal3<S^5 — 'Pa/s'P-yS " '^q.^^^/j + VasV^p) (B13) 

{Sal3<S^5 + Va/sV-^fS — <Sa5<S-y/3 — VasV^p) (B14) 

(TapT^S — TasT^p) . (B15) 



In particular, tree level vertices for these operators are obtained by setting G{p) = I. 

To compute renormalization factors, we project renormalization condition in Eq. |2H1 For 
the axial vector operator, Eq. ISHlis obtained by taking the trace with the tree vertex 757^. 
In the same way, for the four-quark operators, the amputated 4-point Green's function in 
the renormalization condition, Eq. |2Hlis projected into 



A 



jk 



ip) 



a/3;7(5-C//3a;57 



p(4)latt 
j 

-r(4)tree/ \ pfc 

i- j {P)aP;-yS^pa;5^y 



(B16) 

(BIT) 



where the external vertex Ea 



:/3;7<5 

(7m) a/3 (7^)75 + il5li^)ap{l5lfj,)'yS 
{l^l)apil^^)■yS - ilFJl^,)apil5l^^)-yS 
Eap-aS = 5ap5^S - (75) a/3 (75)7,5 

^apS-yS + (75)0/3(75)75 
{<^^iu)ap{<^fiu)^S 

In particular, for the tree level vertex, we obtain 

32/3 

8 -4/3 

AT = 144 X -4/3 2 

5/3 1 

1 7 



(B18) 



(B19) 
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APPENDIX C: RGI FACTOR 



To absorb the momentum dependence of z^^^^*^^ we calculated on the lattice, we use 



' '^scheme 



where the label "scheme" represents the scheme (RI/MOM or MS) and 
asiiJ,) 



(Nf) 2 
MS 



, A lnln(/.VAS^^) 
/5o ln(/i7A.'" 



MS 



(CI) 



(C2) 



7o = 4 

33 - 2Nj 
- ^ 

g 

/3i = 102 - lONf - -Nf. 

3 



Perturbative results for jj-^^n^ were calculated to NLO in Ref. |4£l for both the MS and 



(C3) 
(C4) 

(C5) 



RI/MOM schemes: 



(Nf) 



RI/MOM 



"^MS 



17397- 2070 A^/ + 104iV| 
6(33 - 2Nf)^ 
13095 - 1626Nf + 8N] 
6(33 - 2Nfy 



+ 81n2 



(C6) 
(C7) 
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TABLE I: Basic simulation parameters and numbers of configurations for each observable. Ob- 
servables with the number of configurations used given in bold type are new calculations for this 
paper and the asterisk means the observable was constructed from quark propagators which were 
an average of quark propagators with periodic and anti-periodic boundary conditions in the time 
direction. References refer to previous calculations of the RBC Collaboration. 

DBW2 f3 =1.22 DBW2 (3 =1.04 Wilson f3 =6.0 



Parameters: 

lattice size 24^ x 48 16^ x 32 16^ x 32 

Ls 10 16 16 

Ms 1.65 1.70 1.80 
Configurations used: 

mps 106* 202* 400* 

mv 106* 405 [29j 85 [29] 

fx 106* 202* 400* 

mres, Za 106 405 [29] 85 [29] 

NPR 53 50 40 [21] 

Bk, K-K me 106* 202* 400* [21] 
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TABLE IL Quark masses, nif, used for the measurements hsted in Tabled 



DBW2 /3 = 


1.22 (all) 


0.008, 


0.016, 0.024, 


0.032, 


0.040 






DBW2 (3 = 


1.04 [29] 


0.010, 


0.015, 0.020, 


0.025, 


0.030, 


0.035, 


0.040 




(NPR) 


0.020, 


0.030, 0.040, 


0.050 










(other) 


0.010, 


0.020, 0.030, 


0.040, 


0.050 






Wilson P = 


= 6.0 [29] 


0.010, 


0.015, 0.020, 


0.025, 


0.030, 


0.035, 


0.040 




[21] 


0.010, 


0.020, 0.030, 


0.040, 


0.050 







TABLE IIL Results of the chiral extrapolations of mresi i^^ps fny- We use the linear function 



Co + cinif for m^-cs and cq + Ci(mj + rrires) for rripg and my. The value of rures at nif = 0.02 for 



Wilson 13 = 6.0, and the result of the my = cq + ci mj fit are quoted from Ref. 



2£|. 







Co 


Cl 


xVdof 


mres DBW2 (3 -- 


=1.22 


9.722(27) • 10 


9.5(4.4) • 10-6 


0.637 


DBW2 (3 -- 


=1.04 


1.86(12) • 10" 


5 -4.3(3.0) • 10-5 


0.033 


Wilson /? 


=6.0 


1.24(5) • 10"^ 


(at m/ = 0.02) 




my DBW2 13 -- 


=1.22 


0.2636(48) 


2.418(10) 


0.007 


DBW2 13 -- 


=1.04 


0.3885(59) 


2.34(12) 


0.016 


Wilson (3 


=6.0 


0.404(8) 


2.78(11) 


0.48 


mis DBW2 l3 -- 


=1.22 


0.00142(94) 


1.849(28) 


0.240 


DBW2 (3 -- 


=1.04 


0.0058(15) 


2.584(34) 


1.223 


Wilson P 


=6.0 


0.0052(10) 


3.233(23) 


0.945 
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TABLE IV: Results for my, mps and mps/mv with degenerate and non-degenerate quark masses 
for DBW2 p = 1.22. 



m/ 


my 


mps/mv 


mps 


0.008 


0.2825(55) 


0.462(13) 


0.1304(29) 


0.016 


0.3028(34) 


0.5806(91) 


0.1758(21) 


0.024 


0.3220(27) 


0.6623(75) 


0.2132(18) 


0.032 


0.3412(27) 


0.7217(65) 


0.2462(16) 


0.040 


0.3605(21) 


0.7663(56) 


0.2763(16) 


0.008, 0.016 


0.2928(42) 


0.529(11) 


0.1549(25) 


0.008, 0.024 


0.3029(36) 


0.5836(99) 


0.1768(23) 


0.008, 0.032 


0.3130(33) 


0.6284(93) 


0.1967(22) 


0.008, 0.040 


0.3232(32) 


0.6658(88) 


0.2152(21) 


0.016, 0.024 


0.3125(30) 


0.6251(83) 


0.1953(19) 


0.016, 0.032 


0.3223(28) 


0.6629(77) 


0.2136(19) 


0.016, 0.040 


0.3322(26) 


0.6949(73) 


0.2308(18) 


0.024, 0.032 


0.3316(25) 


0.6943(70) 


0.2302(17) 


0.024, 0.040 


0.3414(24) 


0.7218(66) 


0.2464(17) 


0.032, 0.040 


0.3509(22) 


0.7455(60) 


0.2616(16) 
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TABLE V: Results for my, mps and mps/my quoted from Refs. |2g| and |2g| for DBW2 P = 
1.04 and Wilson (3 = 6.0, my, respectively. Only values in column 6 are newly calculated from 
simultaneous fits, enforcing a common value of mps, to the pseudoscalar-axial and pseudoscalar- 
pseudoscalar correlation functions on the doubled lattice. 





rrif 
J 


my 


mps/my 


™PS (P-P correl.) 


mps (simul. fit) 




0.010 


0.4132(56) 


0.4258(71) 


0.1794(22) 


0.1823(38) 


DBW2 


0.020 


0.4351(37) 


0.5422(50) 


0.2377(15) 


0.2379(29) 


P = 1.04 


0.030 


0.4586(29) 


0.6229(39) 


0.2868(12) 


0.2864(25) 


[29] 


0.040 


0.4825(25) 


0.6825(33) 


0.3300(11) 


0.3298(23) 




0.050 








0.3697(21) 




0.010 


0.442(10) 




0.203(3) 


0.2058(15) 


Wilson 


0.020 


0.462(6) 




0.270(3) 


0.2711(12) 


f3 = 6.0 


0.030 


0.488(5) 




0.324(2) 


0.3245(11) 


m 


0.040 


0.515(4) 




0.371(2) 


0.3716(10) 




0.050 








0.4147(10) 



TABLE VL Results of the fit of to the quenched chiral expansion (Eq.Hg) for DBW2 f3 = 1.22. 
In particular, the results in the final three lines of the table are obtained from a fully covariant fit. 



range 


of mj 


«7r 


6 


b 


xVdof 


0.008 - 


- 0.040 


1.863(27) 


0.028(23) 





1.17 


0.008 - 


- 0.032 


1.804(38) 


0.056(28) 





0.87 


0.016 - 


- 0.032 


1.860(70) 


0.02(4) 





0.13 


0.016 - 


- 0.040 


1.579(96) 


0.093(46) 


5.9(1.2) 


5.8 


0.008 - 


- 0.040 


1.540(98) 


0.101(48) 


6.5(1.2) 


11 
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TABLE VII: Parameters and physical results from meson fits and chiral extrapolations. Results 
for DBW2 (3 = 1.04 are obtained by an extended analysis in Ref. |29j]. Results for Wilson /3 = 6.0 
are quoted from Ref. i29]. 

[GeV] TTT-res [MeV] ams/2 (mK* /mp)iatt J-parameter 
DBW2 (3 = 1.22 2.914(54) 0.2833(54) 0.01474(69) 1.1276(69) 0.387(14) 
DBW2 /3 = 1.04 1.982(30) 0.0368(24) 0.02214(71) 1.138(11) 0.377(37) 
Wilson P = 6.0 1.922(40) 2.38(10) - - - 



TABLE VIII: The square of the pseudoscalar mass, "ipg, expressed in physical units [GeV^]. The 
data for DBW2 (3 = 1.22 reflect jackknife errors while the errors for mpg in lattice units and 
given in Table. I VI II are combined in quadrature for the DBW2, /3 = 1.04 and Wilson /? = 6.0 
results. 





DBW2 (3 = 1.22 


ruf 


DBW2 /3 = 1.04 


Wilson P = 


0.008 


0.1445(80) 


0.01 


0.1305(34) 


0.1565(41) 


0.016 


0.262(11) 


0.02 


0.2224(44) 


0.2715(66) 


0.024 


0.386(15) 


0.03 


0.3221(57) 


0.3889(93) 


0.032 


0.515(20) 


0.04 


0.4272(71) 


0.510(12) 


0.040 


0.648(24) 


0.05 


0.5368(87) 


0.635(15) 
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TABLE IX: Bare values of the pseudoscalar decay constants for each parameter set. Results for 
the three types of analysis discussed in the text, Eqs.ElElandEOl are listed. 





mf 


f(i) 
/ps 


f(2) 

/ps 


f(3) 

/ps 




U.UUo 


U.Uddod( i ( ) 


U.U0i34/(/yj 


U.Uooooi^oD ) 


JJjD WZ 


U.UiO 


n n KCi'7'7 ( '7K\ 
U.Uoy/ /v/oj 


U.UoyzHb/ j 


U.UoUUH j 


P =1.11 


n no A 


U.UdoIU^ /Dj 


U.UOZ r i(DZ j 


n n/^Q /I n/' '7^^ 
U.Uoo4U^ (\) ) 




U.UoZ 




U.UD0y4^Do j 






0.040 


0.06919(83) 


0.06895(68) 


0.06957(85) 




U.Ui 


U.Uo4iU(^yy j 


U.UoziO^yo j 


n CiQQ'Vf 1 oA 

U.Uoo 1 (izj 


L)D VV Z 


n no 
U.Uz 


U.Uoo4o(^oy j 


U.Uo ( 44( (Oj 


U.UooU(^ii ) 


p — i.U4 


n HQ 
U.Uo 


n noo^n^'nQA 
U.Uyzoy^yo j 


u.uyzu (\( Z) 


U.Uyzi^ii j 




U.U4 


u.uyo ( 4^y ( ) 


u.uyoou^ to) 


u.uyoz^ii j 




0.05 


0.1008(10) 


0.10075(76) 


0.1003(12) 




0.01 


0.10043(88) 


0.09920(70) 


0.10006(94) 


Wilson 


0.02 


0.10586(79) 


0.10518(59) 


0.10587(84) 


(3 =6.0 


0.03 


0.11136(79) 


0.11083(58) 


0.11145(84) 




0.04 


0.11672(80) 


0.11628(60) 


0.11682(86) 




0.05 


0.12184(82) 


0.12150(64) 


0.12195(87) 



TABLE X: The renormalization factor Za, linear fitting parameters for the bare pseudoscalar decay 
constants /ps = /o + /i^-ps find bare values of decay constants /^^**^ and f^^^^^\ The statistical 
errors in are taken into account by the jackknife method for DBW2 /? = 1.22 and by quadrature 
for DBW2 p = 1.04 and Wilson p = 6.0. 

ZAjruf = -m,es) /o [GeV] /i xVdof fj'^''^ [GeV] fp''^ [GeV] 

DBW2 /? =1.22 0.88813(19) 0.1547(37) 0.0738(57) 0.033 0.1560(36) 0.1728(37) 
DBW2 (3 =1.04 0.84019(17) 0.1567(33) 0.0813(53) 0.043 0.1582(32) 0.1767(32) 
Wilson P =6.0 0.7555(3) 0.1799(42) 0.0865(38) 0.014 0.1815(42) 0.2011(45) 
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TABLE XI: Results for the physical decay constants obtained from each ensemble and in the 
continuum limit. For the continuum extrapolation, three points with the DBW2 gauge action 
(/? = 1.22, 1.04 and 0.87) are used. 





fn [GeVl 

JUL J 


fK [GeVl 


fx/fn 


DBW2 /3 =1.22 


0.1386(32) 


0.1534(33) 


1.1073(90) 


DBW2 (3 =1.04 


0.1329(27) 


0.1484(27) 


1.1166(85) 


DBW2 (3 =0.87 [29] 


0.1304(67) 


0.1489(52) 


1.142(26) 


Wilson /3 =6.0 


0.1371(32) 


0.1519(34) 


1.1081(50) 


continuum (linear fit) 


0.1395(41) 


0.1528(38) 


1.098(13) 


continuum (constant fit) 


0.1348(20) 


0.1502(19) 


1.114(60) 


(x^/dof for constant fit) 


0.905 


1.128 


0.741 



TABLE XIL The quantity (PS 
ensemble. 



q(A5 2) pg\ -j^ units of [10 ^] for each value of m/ and each 



ruf 


DBW2 /5 = 1.22 


ruf 


DBW2 /3 = 1.04 


Wilson P = 


0.008 


0.708(57) 


0.01 


2.55(19) 


5.57(27) 


0.016 


1.652(84) 


0.02 


5.88(33) 


12.62(48) 


0.024 


2.96(12) 


0.03 


10.48(51) 


21.80(76) 


0.032 


4.60(18) 


0.04 


16.35(72) 


33.2(1.1) 


0.040 


6.59(25) 


0.05 


23.53(98) 


47.0(1.5) 
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TABLE XIII: Fitting parameters for tlie (PS Q(^'^=2) 



PS ) matrix element in lattice units for the 



three types of fitting functions, given by Eas. HUl and HTl 





fit. func. 


ao [10-5] 


ai [10-4GeV-2] 


aa [10-3 GeV-^] 


as 


xVdof 


ag/ai [GeV2] 


DBW2 


(EH) 


0.0 


2.06(23) 


1.047(67) 


-6ai/(47r/)2 


0.073 


-2.013 


P =1.22 


601) 


0.60(53) 


1.83(18) 


1.095(52) 


-6ai/(47r/)2 


0.0001 


-2.013 




m 


0.0 


2.67(71) 


1.035(74) 


-0.25(12) • 10-3 


0.003 


-0.92(68) 


DBW2 


(EH) 


0.0 


8.84(96) 


5.78(31) 


-6ai/(4^/)2 


0.034 


-2.191 


/3 =1.04 


dm) 


-1.5(1.9) 


9.51(94) 


5.59(25) 


-6ai/(4^/)2 


0.010 


-2.191 




m 


0.0 


7.3(2.8) 


5.73(27) 


-2.44(64) • 10-3 


0.018 


-3.3(2.0) 


Wilson 


(El 


0.0 


14.37(82) 


8.08(25) 


-6ai/(47r/)2 


0.058 


-2.057 


(3 =6.0 


(001) 


2.2(2.2) 


13.55(78) 


8.26(18) 


-6ai/(4^/)2 


0.016 


-2.057 




iU) 


0.0 


16.6(2.7) 


8.05(27) 


-2.31(53) • 10-3 


0.007 


-1.39(53) 



TABLE XIV: The pseudoscalar i?-parameter Bps for each ensemble. 





DBW2 13 = 1.22 


ruf 


DBW2 p = 1.04 


Wilson /3 = 6.0 


0.008 


0.499(19) 


0.01 


0.462(18) 


0.505(11) 


0.016 


0.565(11) 


0.02 


0.546(11) 


0.5856(66) 


0.024 


0.6129(81) 


0.03 


0.6016(88) 


0.6361(50) 


0.032 


0.6476(69) 


0.04 


0.6413(71) 


0.6717(41) 


0.040 


0.6740(61) 


0.05 


0.6715(61) 


0.6990(35) 
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TABLE XV: Results for two types of chiral fit to -Bps using the fitting functions in Eqs. l4l-{l and 1441 





fit. func. 


^0 


?i fGeV-^l 


^2 [GeV"2l 


vVdof 


6/^0 fGeV-^l 


n(latt) 
K 


DBW2 
P = 1.22 


m 


0.2719(86) 
0.357(43) 


0.381(20) 
0.344(33) 


-6Co/(4vr/)2 
-0.333(79) 


0.869 
0.006 


-2.013 
-0.93(33) 


0.551(10) 
0.556(12) 


DBW2 
(3 = 1.04 


m 


0.2612(90) 
0.259(39) 


0.407(25) 
0.407(25) 


-6Co/(4vr/)2 
-0.577(86) 


0.033 
0.112 


-2.191 
-2.23(66) 


0.558(10) 
0.558(10) 


Wilson 
/3 = 6.0 


JH 


0.269(14) 
0.293(39) 


0.403(38) 
0.401(40) 


-6Co/(4vr/)2 
-0.482(81) 


0.661 
0.020 


-2.057 
-1.65(47) 


0.5657(63) 
0.5683(70) 



TABLE XVI: Lattice S-parameters for the wrong chirality operators Oj for DBW2 /3 = 1.04. 



m/ 


Byy-AA 


Bss-PP 


Bss+pp 




0.01 


-17.29(54) 


26.62(86) 


-26.9(1.5) 


-13.34(75) 


0.02 


-10.39(19) 


15.54(30) 


-10.27(35) 


-5.13(18) 


0.03 


-7.67(10) 


11.18(16) 


-6.57(16) 


-3.305(79) 


0.04 


-6.193(67) 


8.82(10) 


-4.990(89) 


-2.526(46) 


0.05 


-5.272(48) 


7.347(73) 


-4.111(59) 


-2.094(30) 
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TABLE XVII: The kaon S-parameter Bk and RGI value Bk for each parameter set and choice 
of non-perturbative (NPR) or perturbative (PR) renormahzation factors which are presented in 
Eqs. Kffil and \'A7\ respectively. The continuum values in the first line are obtained from a linear fit 
using the DBW2 p = 1.22 and 1.04 data points. Results for B^ ^^^{fi = 2 GeV) from Wilson 
(3 = 6.0 is reproduced from Ref. |^ for comparison. In the last line, we present the values of 
X^/dof for the constant fit to the dependence. 





NPR 


= 2 GeV) 
PR 


Bk 

NPR 


PR 


continuum (linear fit) 


0.563(21) 


0.547(21) 


0.786(31) 


0.761(28) 


continuum (constant fit) 


0.5357(74) 


0.5342(71) 


0.747(11) 


0.7419(84) 


DBW2 (3 = 1.22 


0.546(11) 


0.539(10) 


0.760(15) 


0.750(14) 


DBW2 (3 = 1.04 


0.526(10) 


0.5298(99) 


0.731(16) 


0.738(10) 


Wilson p = 6.0 


0.535(6) 








(x^/dof for constant fit) 


1.785 


0.404 


1.899 


0.517 
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TABLE XVIII: Estimates of the size of various errors that are present and corrections that are 
required by our result for S^^^^^. This estimate of the finite volume correction, (a), was taken 
from the work reported in Ref. jsij] and the estimate of that for the use of degenerate valence 
quarks quenched approximation, (b), from Ref. Our quoted estimate of all systematic errors 
but quenching (±0.038) is determined by combining the numbers in the second through tenth rows 
in quadrature. 



Source of error 


Estimated size 


(%) 


Statistical 


±0.021 


(4%) 


Excited state contamination 


±0.017 


(3%) 


Determination of fx 


±0.005 


(1%) 


Mass of K meson 


±0.042 


(3%) 


Non-perturbative renormalization 


±0.005 


(1%) 


Perturbative renormalization 


±0.023 


(4%) 


Mixing with wrong chirality operators 


±0.0001 


(0.02%) 


Finite volume 


±0.01" 


(2%) 


Degenerate valence quarks 


±0.015'' 


(3%) 


Continuum extrapolation 


±0.001 


(0.2%) 


Quenching 


±0.030 


(6%) 
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Wilson I 5k sweeps 



DBW2 



OA o- 



FIG. 1: A schematic of our algorithm for generating gauge configurations with the DBW2 action 
and P = 1.22. We first produce gauge configurations using the Wilson gauge action with /? = 6.25 
(vertical direction) and save them every 10,000 sweeps (open circles). These saved configurations 
are used as the initial configuration for a subsequent DBW2 evolution with /3 = 1.22 (horizontal 
direction) . 
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o 



o 



o 




10 20 30 40 
#sweep (x1 Ok) 



10 15 



FIG. 2: The time history of the topological charge for the vertical direction in Fig. ^ The top 
panel is for the initial gauge configurations with Wilson gauge action and (3 = 6.25, given by the 
open circles in Fig. ^ The middle panel is for the DBW2 action after 5,000 heatbath sweeps and 
the bottom is for the DBW2 case after 10,000 sweeps. For each panel, the histogram in the right 
of the panel shows the distribution of Qtop- 
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3.20e-03 

3.10e-03 
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• m,=0.00 
O m,=0.02 
□ m,=0.03 
O m,=0.04 




01 23456789 
#sweep [1 k] 



FIG. function of the number of sweeps for the DBW2 evolutions with (3 — 1.22. mres 

was measured for = 8 and M5 = 1.7 and these results are averaged over 20 of the 53 DBW2 
evolutions. 




FIG. 4: The pseudoscalar meson mass, mps, as a function of the number of sweeps in the DBW2 
evolution. 
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5 10 15 20 25 

t 

FIG. 5: Ratio of the right-hand side of Eq. IHlfor the DBW2 (3 = 1.22 data as a function of t. From 
bottom to top panel, data are plotted in order of hghtness of the quark. 
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FIG. 6: Residual quark mass, rrij-es, for DBW2 /? = 1.22. The value at mf = is obtained from a 
linear extrapolation. 
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FIG. 7: The pscudoscalar effective mass for the DBW2, [3 = 1.22 data set plotted as a function of 
the time t. This was obtained from the point-wah correlator C^^. 
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FIG. 8: Vector meson effective mass plots for the DBW2, (3 = 1.22 ensemble. Lines denote fitting 
results and indicate fitting range, central values and jackknife errors. 




FIG. 10: Pseudoscalar meson mass-squared for DBW2 (3 = 1.22. Values computed from non- 
degenerate quark masses mi ^ m-2 are plotted as mf + m^es = {n^i + 1^2 + 2mres)/2. 
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FIG. 11: The J-parameter for DBW2 /3 = 1.22 (filled circle), 1.04 (filled square) and 0.87 (filled 
triangle). The latter is reproduced from Ref. 29]. The open triangle represents the value obtained 
with Wilson f3 = 6.0. Scaling violations are evidently less than the statistical errors. 
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FIG. 12: Deviation of /^g* from /j,^] as a function of rn-pg [GeV^] for each parameter set, DBW2 
(3 = 1.22 (filled circles), /3 = 1.04 (filled squares), and Wilson f3 = 6.0 (open circles). 
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FIG. 13: Effective lattice decay constant fp^ obtained from the pseudoscalar effective mass and 
corresponding amplitudes for DBW2 (3 = 1.22. Tfie mass ruf increases from bottom to top. The 
fitting range 18 < t < 31 was used to extract fpg for this case. 
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FIG. 14: Same plot as Fig.UHlbut for DBW2 /? = 1.04. The fitting range 12 < t < 19 was used to 
extract fp^ for this case. 
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FIG. 15: Renormalization factor Za as a function of mj + rrires for DBW2 /5 = 1.22. 



68 



0.20 
0.19 
0.18 
0.17 

> 

a 0.16 

in 
a. 

0.15 
0.14 
0.13 
0.12 

0.1 0.2 0.3 0.4 0.5 0.6 0.7 
mp/ [GeV^] 

FIG. 16: Renormalized values of the pseudoscalar decay constant as a function of rrip^ [GeV^] for 
DBW2 /? = 1.22 (filled circle), /3 = 1.04 (filled square), and Wilson /3 = 6.0 (open diamond). 
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FIG. 17: Scaling property of the ratio of physical pseudoscalar decay constants, fK/fn- The values 
for P = 0.87 (filled diamond) and /? = 1.04 (open rectangle) are taken from Ref. 29]. The linear 
extrapolation to the continuum limit (solid line) uses the three new DBW2 data points: (3 = 1.22 
(filled circle), f3 = 1.04 (filled square) and /3 = 0.87 (filled diamond). 
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FIG. 18: Scaling properties of the decay constants and fx- The values for /? = 0.87 (triangles) 
are taken from Ref. 3]. Linear extrapolations to the continuum limit (solid lines) use the three 
DBW2 data points, (3 = 1.22 (circles), (3 = 1.04 (squares) and /3 = 0.87 (triangles). Experimental 
results are also shown. 



0.92 
0.90 
0.88 
0.86 
0.84 
0.82 
0.80 
0.78 
0.76 




0.0 



0.5 



• m,=-m,^j 
O m,=0.008 
□ m,=0.016 
O m,=0.024 
A m,=0.032 
<! m,=0.040 



1.0 



1.5 



2.0 



2.5 
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the open symbols are linearly extrapolated to the value (filled circles) in the chiral limit. 



70 



0.01 



0.00 



-0.01 



-0.02 



-0.03 



VV+AA; VV-AA,p=1.22 





• m,=-m^|,j 




O m,=0.008 




□ m,=0.016 




O m,=0.024 




A m,=0.032 


1 

1 , 


< m,=0.040 



0.0 



0.5 



1.0 1.5 



2.0 



2.5 



0.01 



0.00 



-0.01 



-0.02 



-0.03 




FIG. 20: Same figure as Fig. but for the renormalization factor Z "^Zvv+AA, vv-AA 
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FIG. 21: Renormalization factors "^Zyv+AA, vv-AA as a function of -l-mres for several fixed 
values of p^jj^.^, for each of which a linear fit is drawn. 
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FIG. 22: Renormalization factors '^Zyy+AA, j as a function of -t-mres for several fixed values 
of pf^tf From top to bottom, j = SS - PP, SS + PP and TT . Results for DBW2 /3 = 1.22 and 
1.04 are in the left and right column, respectively. 
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FIG. 23: Renormalization factor of Bk as a function of pf^^^^ for DBW2 /3 = 1.22 (left) and (3 = 1.04 
(right). In each panel, Zbj^ (circles), w~^{pu\.t)ZBK (squares) and its linear extrapolation using 
data for piatt > 1 are shown. 
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FIG. 24: The time-dependence of the numerator of the left-hand-side of Eq. |2H1 for the DBW2 
j3 = 1.22 data set. The plotted horizontal lines indicate both the fitting range and the upper and 
lower limits of the resulting fitting value. 
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FIG. 25: The same quantity as plotted in Fig. 01 except for the DBW2 /3 = 1.04 data set. 
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FIG. 27: Time dependence oi Bps for /? = 1.22. Panels from bottom to top correspond to increasing 
ruf. The horizontal bars in each plot indicate the range used in the constant fit, the results for the 
central value (solid line) and the jackknife error (dashed lines). 
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FIG. 28: Time dependence of Bps for /3 = 1.04. The organization is same as in Fig. 1271 



75 



0.7 



0.6 - 



• DBW2 p=1.22 

O mps=m^ 

— constrainted chral log. 
- free chiral log. 




0.7 



0.6 - 



^ 0.5 



0.4 - 



0.3 



• DBW2 p=1.04 
— constrainted chiral log. 



0.1 0.2 0.3 0.4 0.5 0.6 0.7 

mp/ [GeV'] 



FIG. 29: Bps as a function of mj^g [GeV^] for DBW2 /? = 1.22 (left) and (3 = 1.04 (right). In each 
panel, solid curves are the results of the fit to Eq. Dashed curve in the left panel denotes the 
fit to Eq. m 
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FIG. 30: Bps renormalized in MS NDR scheme at /i = 2 GeV as a function of [GeV^]. Filled 
symbols are from DBW2 /3 = 1.22 (circle) and 1.04 (square). Fitting curves indicated by the solid 
and dashed lines respectively are added to them. Open diamonds are Wilson (3 = 6.0 data. 



76 



♦ 


1 




1 


1 


• VV+AA 
■ VV-AA 


- 










♦ SS-PP 

ASS+PP 

^ T — r 


- 










A. 


-•- 


• 




• 


• 


-•- 




< 




< 

A 
■ 










































1 








1 


1 


0.01 


0.02 




0.03 


0.04 


0.05 



m, 

FIG. 31: B-parameters for Oyv+AA (circle), Ovv-AA (square), Oss-PP (diamond), Oss+pp 
(triangle) and Ott (left-triangle) for DBW2 /3 = 1.04 as a function of m/. 
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FIG. 32: Summary of our results for ^^^(/i = 2 GeV) renormalized with Nf = as a function 
of the lattice spacing squared. The filled circles are our results and the open symbols are quoted 

A'^^ = 2 result obtained in Ref. |3ll |. 
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FIG. 33: The continuum extrapolation of Bps evaluated in the chiral limit. 



